Parte 1: Optimizacion Numerica.¶
pip install deap
Requirement already satisfied: deap in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (1.3.3)Note: you may need to restart the kernel to use updated packages. Requirement already satisfied: numpy in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from deap) (1.22.3)
pip install pyswarm
Requirement already satisfied: pyswarm in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (0.6) Requirement already satisfied: numpy in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pyswarm) (1.22.3) Note: you may need to restart the kernel to use updated packages.
pip install geopandas
Requirement already satisfied: geopandas in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (0.12.2) Requirement already satisfied: shapely>=1.7 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from geopandas) (1.8.0) Requirement already satisfied: pandas>=1.0.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from geopandas) (1.4.2) Requirement already satisfied: fiona>=1.8 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from geopandas) (1.8.20) Requirement already satisfied: pyproj>=2.6.1.post1 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from geopandas) (3.1.0) Requirement already satisfied: packaging in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from geopandas) (21.3) Requirement already satisfied: attrs>=17 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (21.4.0) Requirement already satisfied: certifi in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (2022.12.7) Requirement already satisfied: click>=4.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (7.1.2) Requirement already satisfied: cligj>=0.5 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (0.7.2) Requirement already satisfied: click-plugins>=1.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (1.1.1) Requirement already satisfied: six>=1.7 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (1.16.0) Requirement already satisfied: munch in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (2.5.0) Requirement already satisfied: setuptools in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (61.2.0) Requirement already satisfied: python-dateutil>=2.8.1 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pandas>=1.0.0->geopandas) (2.8.2) Requirement already satisfied: numpy>=1.18.5 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pandas>=1.0.0->geopandas) (1.22.3) Requirement already satisfied: pytz>=2020.1 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pandas>=1.0.0->geopandas) (2022.1) Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from packaging->geopandas) (3.0.9) Note: you may need to restart the kernel to use updated packages.
pip install pygad
Requirement already satisfied: pygad in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (2.19.2) Requirement already satisfied: numpy in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pygad) (1.22.3) Requirement already satisfied: matplotlib in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pygad) (3.5.2) Requirement already satisfied: cloudpickle in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pygad) (2.0.0) Requirement already satisfied: fonttools>=4.22.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (4.33.3) Requirement already satisfied: pillow>=6.2.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (9.0.1) Requirement already satisfied: packaging>=20.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (21.3) Requirement already satisfied: python-dateutil>=2.7 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (2.8.2) Requirement already satisfied: kiwisolver>=1.0.1 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (1.4.3) Requirement already satisfied: pyparsing>=2.2.1 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (3.0.9) Requirement already satisfied: cycler>=0.10 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (0.11.0) Requirement already satisfied: six>=1.5 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from python-dateutil>=2.7->matplotlib->pygad) (1.16.0) Note: you may need to restart the kernel to use updated packages.
#Importamos las librerias Necesarias
import random
import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go
from matplotlib.colors import LogNorm
import matplotlib.animation as animation
from IPython.display import HTML
from deap import algorithms, base, creator, tools
from pyswarm import pso
from scipy.optimize import minimize
import geopandas as gpd
from shapely.geometry import Point
import pyproj
import pygad
from pygad import gann
from matplotlib import animation, rc
import pandas as pd
import itertools
import imageio
import os
from matplotlib.animation import FuncAnimation
# define la función de Rastrigin de dos dimensiones
def f_rastrigin2d(x1, x2):
y = 20 + x1 ** 2 - 10 * np.cos(2 * np.pi * x1) + x2 ** 2 - 10 * np.cos(2 * np.pi * x2)
return y
# define la función de Rastrigin de tres dimensiones
def f_rastrigin3d(x1,x2,x3):
y = 20 + x1 ** 2 - 10 * np.cos(2 * np.pi * x1) + x2 ** 2 - 10 * np.cos(2 * np.pi * x2) + x3 ** 2 - 10 * np.cos(2 * np.pi * x3)
return y
n = 50 # define el número de puntos en la malla
# crea una lista de valores igualmente espaciados desde -5.12 a 5.12
x1 = np.linspace(-5.12, 5.12, n)
x2 = np.linspace(-5.12, 5.12, n)
# crea una malla de valores de entrada 2D utilizando los valores de x1 y x2
X1, X2 = np.meshgrid(x1, x2)
# aplana la malla y concatena en una matriz de entrada X
X = np.vstack((X1.flatten(), X2.flatten())).T
# calcula el valor de la función de Rastrigin de dos dimensiones para cada punto en la matriz X
z = f_rastrigin2d(X[:, 0], X[:, 1])
# reformatea el resultado en una matriz 2D de tamaño n x n
Z = z.reshape(n, n)
# Crear un gráfico de contorno con valores numéricos
levels = np.arange(0, 90, 10)
cp = plt.contourf(x1, x2, Z, levels=levels, cmap='plasma')
plt.colorbar(cp)
# Agregar líneas de contorno
plt.contour(x1, x2, Z, levels=levels, colors='k', linewidths=0.5)
# Título y etiquetas de los ejes
plt.title("Función de Rastrigin")
plt.xlabel("x1")
plt.ylabel("x2")
# Mostrar el gráfico
plt.show()
# Cálculo de la función
z = f_rastrigin2d(X1, X2)
# Creación de la figura
fig = go.Figure(data=[go.Surface(x=X1, y=X2, z=z)])
# Personalización de la figura
fig.update_layout(
title='Función de Rastrigin 2D',
scene = dict(
xaxis_title='X1',
yaxis_title='X2',
zaxis_title='Y',
aspectmode='cube',
camera=dict(
eye=dict(x=1.75, y=-1.75, z=1.25)
)
)
)
# Mostrar la figura
fig.show()
Vamos a usar el metodo de descenso de gradiente llamado metodo de Newton- Raphson
# Define una función que calcula la función de Rastrigin dada una entrada x
def rastrigin(x):
n = len(x) # Obtiene la longitud de la entrada x
return 10*n + np.sum(x**2 - 10*np.cos(2*np.pi*x)) # Calcula y devuelve el valor de la función de Rastrigin
# Define una función que calcula el gradiente de la función de Rastrigin dada una entrada x
def gradient(x):
n = len(x) # Obtiene la longitud de la entrada x
grad = np.zeros(n) # Crea un vector gradiente de ceros de longitud n
for i in range(n): # Itera sobre cada entrada de la entrada x
grad[i] = 2*x[i] + 20*np.pi*np.sin(2*np.pi*x[i]) # Calcula la entrada correspondiente en el vector gradiente
return grad # Devuelve el vector gradiente
# Define una función que calcula la matriz Hessiana de la función de Rastrigin dada una entrada x
def hessian(x):
n = len(x) # Obtiene la longitud de la entrada x
hess = np.zeros((n, n)) # Crea una matriz Hessiana de ceros de tamaño n x n
for i in range(n): # Itera sobre cada entrada de la entrada x
for j in range(n): # Itera sobre cada entrada de la entrada x
if i == j: # Si i y j son iguales
hess[i][i] = 2 + 40*np.pi**2*np.cos(2*np.pi*x[i]) # Calcula la entrada correspondiente en la matriz Hessiana
else:
hess[i][j] = 0 # De lo contrario, establece la entrada correspondiente en la matriz Hessiana a cero
return hess # Devuelve la matriz Hessiana
# Define una función que optimiza la función de Rastrigin dada una conjetura inicial, un número máximo de iteraciones y una tolerancia
def optimize_rastrigin(initial_guess, max_iterations, tolerance=1e-6):
x = initial_guess # Establece la conjetura inicial como el valor actual de x
path = [x] # Crea una lista que almacena los valores de x en cada iteración
for i in range(max_iterations): # Itera sobre el número máximo de iteraciones
grad = gradient(x) # Calcula el gradiente de la función de Rastrigin en el valor actual de x
hess = hessian(x) # Calcula la matriz Hessiana de la función de Rastrigin en el valor actual de x
step = np.linalg.solve(hess, -grad) # Calcula el paso de optimización utilizando la solución de la ecuación Hessiana
x_new = x + step # Calcula el nuevo valor de x utilizando el paso de optimización
if np.linalg.norm(x_new - x) < tolerance: # Si la norma de la diferencia entre el nuevo y el valor actual de x es menor que la tolerancia, la optimización se ha completado
break
x = x_new # De lo contrario, establece el nuevo valor de x como el valor actual de x y continúa iterando
path.append(x) # Agrega el nuevo valor de x a la lista de valores de x en cada iteración
return x, path # Devuelve el valor final
rastrigin_init_point_2D = np.random.uniform(-5.12, 5.12, size=2) #Se genera un número aleatorio y se acota entre el dominio de la función
max_iterations = 1000
rastrigin_min_2D,rastrigin_path_2D = optimize_rastrigin(rastrigin_init_point_2D, 1000)
#Impresión de resultados
print(f"Init value: {rastrigin_init_point_2D}")
print(f"Minimum: {rastrigin_min_2D}")
print(f"Function value at minimum: {rastrigin(rastrigin_min_2D)}")
print("Optimization path: ")
for point in rastrigin_path_2D:
print(f"Punto evaluado:{point}--", f"Valor de la función en el punto: {rastrigin(point)}")
# Generate a grid of x, y values
rtg_x_min, rtg_x_max, rtg_x_maxvalues = -5.12, 5.12, 100
rtg_y_min, rtg_y_max, rtg_y_maxvalues = -5.12, 5.12, 100
x_rastrigin = np.linspace(rtg_x_min, rtg_x_max, rtg_x_maxvalues)
y_rastrigin = np.linspace(rtg_y_min, rtg_y_max, rtg_y_maxvalues)
X_rastrigin, Y_rastrigin = np.meshgrid(x_rastrigin, y_rastrigin)
# Compute the Rastrigin function for each (x, y) pair
Z_rastrigin = np.zeros_like(X_rastrigin)
for i in range(X_rastrigin.shape[0]):
for j in range(X_rastrigin.shape[1]):
Z_rastrigin[i, j] = rastrigin(np.array([X_rastrigin[i, j], Y_rastrigin[i, j]]))
rastrigin_minimum= rastrigin_min_2D.reshape(-1, 1)
#Test plot
rastrigin_surface_fig = plt.figure(figsize=(10,8))
rastrigin_surface_axes = plt.axes(projection='3d', elev=50, azim=-50)
rastrigin_surface_axes.plot_surface(X_rastrigin, Y_rastrigin, Z_rastrigin, norm=LogNorm(), rstride=1, cstride=1,
edgecolor='none', alpha=.2, cmap=plt.cm.jet)
rastrigin_surface_axes.plot(*rastrigin_minimum, rastrigin(rastrigin_minimum), 'b*', markersize=10)
rastrigin_surface_axes.set_xlabel('$x$')
rastrigin_surface_axes.set_ylabel('$y$')
rastrigin_surface_axes.set_zlabel('$z$')
rastrigin_surface_axes.set_xlim((rtg_x_min, rtg_x_max))
rastrigin_surface_axes.set_ylim((rtg_y_min, rtg_y_max))
plt.show()
# Create a 2D contour plot of the Rastrigin function and plot the ptimization path
rastrigin_contour_fig, rastrigin_contour_axes = plt.subplots(figsize=(10,8))
contour = rastrigin_contour_axes.contour(X_rastrigin, Y_rastrigin, Z_rastrigin, levels=20, cmap='coolwarm')
rastrigin_contour_axes.set_xlabel('x')
rastrigin_contour_axes.set_ylabel('y')
for point in rastrigin_path_2D:
rastrigin_contour_axes.plot(*point, 'o-',color='blue',alpha=1)
# Add a marker for the global minimum at (0, 0)
rastrigin_contour_axes.scatter(*rastrigin_minimum, marker='*', s=500, color='red')
plt.colorbar(contour)
plt.show()
Init value: [-0.28367766 -1.40353059] Minimum: [-0.99495911 -1.50764073] Function value at minimum: 23.256417941943447 Optimization path: Punto evaluado:[-0.28367766 -1.40353059]-- Valor de la función en el punto: 32.369204323421876 Punto evaluado:[-1.04988166 -1.5232434 ]-- Valor de la función en el punto: 23.803209023269208 Punto evaluado:[-0.99301162 -1.50755286]-- Valor de la función en el punto: 23.25716802783503 Punto evaluado:[-0.99495911 -1.50764073]-- Valor de la función en el punto: 23.256417941943447
# Define the function to update the animation
def update(frame):
x = rastrigin_path_2D[frame]
line.set_data(x[0], x[1])
if frame > 0:
line_trace.set_data([rastrigin_path_2D[i][0] for i in range(frame+1)], [rastrigin_path_2D[i][1] for i in range(frame+1)])
return line, line_trace
# Plot the contour and the starting point
rastrigin_contour_fig, rastrigin_contour_axes = plt.subplots(figsize=(10,8))
rastrigin_contour_axes.set_title('Optimization Path Animation')
rastrigin_contour_axes.set_xlabel('x')
rastrigin_contour_axes.set_ylabel('y')
rastrigin_contour = rastrigin_contour_axes.contour(X_rastrigin, Y_rastrigin, Z_rastrigin, levels=20, cmap='coolwarm')
rastrigin_contour_axes.plot(rastrigin_init_point_2D[0], rastrigin_init_point_2D[1], 'o', color='blue')
# Set up the animation
line, = rastrigin_contour_axes.plot([], [], 'o-', color='red', lw=10)
line_trace, = rastrigin_contour_axes.plot([], [], '-', color='green', lw=4)
ani = animation.FuncAnimation(rastrigin_contour_fig, update, frames=len(rastrigin_path_2D), interval=500, blit=True)
HTML(ani.to_jshtml())
rastrigin_init_point_3D = np.random.uniform(-5.12, 5.12, size=3) #Se genera un número aleatorio y se acota entre el dominio de la función
rastrigin_min_3D,rastrigin_path_3D = optimize_rastrigin(rastrigin_init_point_3D, 1000)
#Impresión de resultados
print("Init value", rastrigin_init_point_3D)
print("Minimum:", rastrigin_min_3D)
print("Function value at minimum:", rastrigin(rastrigin_min_3D))
print("Optimization path: ")
for point in rastrigin_path_3D:
print(f"Punto evaluado:{point}--", f"Valor de la función en el punto: {rastrigin(point)}")
Init value [-2.87663391 2.66267468 0.51743905] Minimum: [-2.9848557 2.51274332 0.50254604] Function value at minimum: 55.487715409989356 Optimization path: Punto evaluado:[-2.87663391 2.66267468 0.51743905]-- Valor de la función en el punto: 53.64496918876429 Punto evaluado:[-3.01119854 2.42587958 0.50249088]-- Valor de la función en el punto: 54.16324400210692 Punto evaluado:[-2.98482268 2.52015726 0.50254604]-- Valor de la función en el punto: 55.47697082306979 Punto evaluado:[-2.9848557 2.51272392 0.50254604]-- Valor de la función en el punto: 55.487715336295715 Punto evaluado:[-2.9848557 2.51274332 0.50254604]-- Valor de la función en el punto: 55.487715409989356
# Definir los límites de los parámetros
limite_inferior = -5.0
limite_superior = 5.0
# Definir los parámetros del algoritmo genético
num_generaciones = 50
tamano_poblacion = 100
probabilidad_mutacion = 0.2
probabilidad_cruce = 0.5
# Crear los tipos de objeto para el problema
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individuo", list, fitness=creator.FitnessMin)
# Crear la caja de herramientas del algoritmo genético
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, limite_inferior, limite_superior)
toolbox.register("individuo", tools.initCycle, creator.Individuo,
(toolbox.attr_float, toolbox.attr_float), n=1)
toolbox.register("poblacion", tools.initRepeat, list, toolbox.individuo)
# Definir la función de evaluación
def evaluar(individuo):
return f_rastrigin2d(individuo[0], individuo[1]),
# Definir los operadores genéticos
toolbox.register("evaluate", evaluar)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)
# Crear la población inicial
poblacion = toolbox.poblacion(n=tamano_poblacion)
# Ejecutar el algoritmo genético
for generacion in range(num_generaciones):
# Seleccionar a los padres
padres = toolbox.select(poblacion, len(poblacion))
# Clonar los padres para aplicar operadores genéticos
hijos = list(map(toolbox.clone, padres))
# Aplicar cruce y mutación a los hijos
for hijo in hijos:
if random.random() < probabilidad_cruce:
toolbox.mate(hijo, random.choice(hijos))
del hijo.fitness.values
if random.random() < probabilidad_mutacion:
toolbox.mutate(hijo)
del hijo.fitness.values
# Evaluar la aptitud de los nuevos individuos
individuos_nuevos = [individuo for individuo in hijos if not individuo.fitness.valid]
aptitudes_nuevas = toolbox.map(toolbox.evaluate, individuos_nuevos)
for individuo, aptitud in zip(individuos_nuevos, aptitudes_nuevas):
individuo.fitness.values = aptitud
# Reemplazar la población anterior por la nueva generación
poblacion = hijos
# Seleccionar al mejor individuo
mejor_individuo = tools.selBest(poblacion, 1)[0]
print("Mejor solución encontrada:", mejor_individuo)
print("Valor de la función objetivo:", mejor_individuo.fitness.values[0])
Mejor solución encontrada: [8.604753801223136e-09, 3.210492889238655e-09] Valor de la función objetivo: 1.5987211554602254e-14
# Definir los límites de los parámetros
limite_inferior = -5.0
limite_superior = 5.0
# Definir los parámetros del algoritmo genético
num_generaciones = 50
tamano_poblacion = 100
probabilidad_mutacion = 0.2
probabilidad_cruce = 0.5
# Crear los tipos de objeto para el problema
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individuo", list, fitness=creator.FitnessMin)
# Crear la caja de herramientas del algoritmo genético
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, limite_inferior, limite_superior)
toolbox.register("individuo", tools.initRepeat, creator.Individuo,
toolbox.attr_float, n=3)
toolbox.register("poblacion", tools.initRepeat, list, toolbox.individuo)
# Definir la función de evaluación
def evaluar(individuo):
return f_rastrigin3d(individuo[0], individuo[1], individuo[2]),
# Definir los operadores genéticos
toolbox.register("evaluate", evaluar)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)
# Crear la población inicial
poblacion = toolbox.poblacion(n=tamano_poblacion)
# Ejecutar el algoritmo genético
for generacion in range(num_generaciones):
# Seleccionar a los padres
padres = toolbox.select(poblacion, len(poblacion))
# Clonar los padres para aplicar operadores genéticos
hijos = list(map(toolbox.clone, padres))
# Aplicar cruce y mutación a los hijos
for hijo in hijos:
if random.random() < probabilidad_cruce:
toolbox.mate(hijo, random.choice(hijos))
del hijo.fitness.values
if random.random() < probabilidad_mutacion:
toolbox.mutate(hijo)
del hijo.fitness.values
# Evaluar la aptitud de los nuevos individuos
individuos_nuevos = [individuo for individuo in hijos if not individuo.fitness.valid]
aptitudes_nuevas = toolbox.map(toolbox.evaluate, individuos_nuevos)
for individuo, aptitud in zip(individuos_nuevos, aptitudes_nuevas):
individuo.fitness.values = aptitud
# Reemplazar la población anterior por la nueva generación
poblacion = hijos
# Seleccionar al mejor individuo
mejor_individuo = tools.selBest(poblacion, 1)[0]
print("Mejor solución encontrada:", mejor_individuo)
print("Valor de la función objetivo:", mejor_individuo.fitness.values[0])
Mejor solución encontrada: [0.019417809579997895, 0.009269224741063526, 0.0009761931431016592] Valor de la función objetivo: -9.908058456113144
c:\Users\Hp\anaconda3\envs\tf_gpu\lib\site-packages\deap\creator.py:138: RuntimeWarning: A class named 'FitnessMin' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it. c:\Users\Hp\anaconda3\envs\tf_gpu\lib\site-packages\deap\creator.py:138: RuntimeWarning: A class named 'Individuo' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it.
#Definimos la funcion rastringin para los calculos
def f(x, y):
return (20 + x ** 2 - 10 * np.cos(2 * np.pi * x) + y ** 2 - 10 * np.cos(2 * np.pi * y))
#Calculamos los datos X Y de la funcion
x, y = np.array(np.meshgrid(np.linspace(-5.12, 5.12, 50),
np.linspace(-5.12, 5.12, 50)))
#Evalamos los datos X Y
z = f(x, y)
#Encontramos el minimo Global de la funcion
x_min = x.ravel()[z.argmin()]
y_min = y.ravel()[z.argmin()]
# Se inicializan los Hiperparametros del algoritmo PSO
c1 = c2 = 0.1 # Factores de Acelearacion
w = 0.8 # Coeficiente de Inercia
# Creamos las particulas
n_particles = 35
np.random.seed(100)
X = np.random.rand(2, n_particles) * 5 # Posiscion Inicial de las particulas
V = np.random.randn(2, n_particles) * 0.5 # velocidad de las particulas
# Inicializacion de datos
pbest = X # Mejor Posicion de la particula
pbest_obj = f(X[0], X[1])
gbest = pbest[:, pbest_obj.argmin()] # Mejor Posicion del Enjambre
gbest_obj = pbest_obj.min()
#funciond de actualizacion de los puntos del enjambre
def update():
global V, X, pbest, pbest_obj, gbest, gbest_obj
# actualizacion de parametros
r1, r2 = np.random.rand(2)
# Actualizacion de velocidades
V = w * V + c1*r1*(pbest - X) + c2*r2*(gbest.reshape(-1, 1)-X)
X = X + V # Actualizacion de Posicion de particulas
obj = f(X[0], X[1]) # Nueva Posicion
# Se evalua si la posicion nueva 'obj' es mejor a la anterior
pbest[:, (pbest_obj >= obj)] = X[:, (pbest_obj >= obj)]
pbest_obj = np.array([pbest_obj, obj]).min(axis=0)
gbest = pbest[:, pbest_obj.argmin()]
gbest_obj = pbest_obj.min()
# Visualizacion del grafico de Grienwank Y los puntos del enjambre
fig, ax = plt.subplots(figsize=(16, 8))
fig.set_tight_layout(True)
img = ax.imshow(z, extent=[-5.12, 5.12, -5.12, 5.12],
origin='lower', cmap='plasma', alpha=0.5)
fig.colorbar(img, ax=ax)
ax.plot([x_min], [y_min], marker='x', markersize=10,
color="white") # Minimo Global
contours = ax.contour(x, y, z, 10, colors='black', alpha=0.4)
ax.clabel(contours, inline=True, fontsize=8, fmt="%.0f")
# mejor posision del punto
pbest_plot = ax.scatter(
pbest[0], pbest[1], marker='o', color='black', alpha=0.5)
p_plot = ax.scatter(X[0], X[1], marker='o', color='green', alpha=0.5) # Puntos
p_arrow = ax.quiver(X[0], X[1], V[0], V[1], color='green', width=0.005,
angles='xy', scale_units='xy', scale=1) # flecha de los puntos
gbest_plot = plt.scatter([gbest[0]], [gbest[1]], marker='*',
s=100, color='red', alpha=0.4) # Mejor posicion del enjambre
ax.set_xlim([-5.12, 5.12])
ax.set_ylim([-5.12, 5.12])
def animate(i):
# Animacion de los puntos en moviento
title = 'Iteration {:02d}'.format(i)
# Update params
update()
# Set picture
ax.set_title(title)
pbest_plot.set_offsets(pbest.T)
p_plot.set_offsets(X.T)
p_arrow.set_offsets(X.T)
p_arrow.set_UVC(V[0], V[1])
gbest_plot.set_offsets(gbest.reshape(1, -1))
return ax, pbest_plot, p_plot, p_arrow, gbest_plot
anim = FuncAnimation(fig, animate, frames=list(
range(1, 50)), interval=500, blit=False, repeat=True)
# Se guarda la animacion en formato gif
anim.save("PSORAS.gif", dpi=120, writer="imagemagick")
# Informacion genreal del algoritmo luego de optimizar
print("PSO found best solution at f({})={}".format(gbest, gbest_obj))
print("Global optimal at f({})={}".format([x_min,y_min], f(x_min,y_min)))
MovieWriter imagemagick unavailable; using Pillow instead.
PSO found best solution at f([-0.0005921 -0.00046606])=0.00011264565412716365 Global optimal at f([-0.9404081632653059, -0.9404081632653059])=3.1543848912857424
from IPython.display import Image, display
Image(filename='PSORAS.gif')
<IPython.core.display.Image object>
lb = [-10, -10, -10] # límite inferior de las variables
ub = [10, 10, 10] # límite superior de las variables
x_opt, f_opt = pso(lambda x: f_rastrigin3d(x[0], x[1], x[2]), lb, ub)
print("El punto óptimo es:", x_opt)
Stopping search: Swarm best objective change less than 1e-08 El punto óptimo es: [ 2.05480392e-05 -9.57955878e-06 1.71508643e-05]
# Definir los parámetros de la evolución diferencial
F = 0.5 # factor de amplificación
CR = 0.9 # tasa de recombinación
NP = 20 # tamaño de la población
ngen = 100 # número de generaciones
# Definir el límite inferior y superior para los valores de x e y
lim_inf = -10
lim_sup = 10
# Generar la población inicial de forma aleatoria
pop = []
for i in range(NP):
x = random.uniform(lim_inf, lim_sup)
y = random.uniform(lim_inf, lim_sup)
pop.append((x, y))
# Realizar la evolución diferencial
for gen in range(ngen):
for i in range(NP):
# Seleccionar 3 individuos de forma aleatoria
a, b, c = random.sample(pop, 3)
# Seleccionar un número aleatorio para cada variable
R1, R2 = random.random(), random.random()
# Realizar la mutación diferencial
x_new = np.clip(a[0] + F*(b[0] - c[0]), lim_inf, lim_sup)
y_new = np.clip(a[1] + F*(b[1] - c[1]), lim_inf, lim_sup)
# Realizar la recombinación
if random.random() < CR:
x = x_new
y = y_new
else:
x = a[0]
y = a[1]
# Evaluar la función en el nuevo individuo
fitness_new = f_rastrigin2d(x, y)
# Reemplazar el individuo anterior si el nuevo individuo es mejor
if fitness_new < f_rastrigin2d(pop[i][0], pop[i][1]):
pop[i] = (x, y)
# Imprimir el mejor individuo encontrado
best_ind = min(pop, key=lambda x: f_rastrigin2d(x[0], x[1]))
print('Valor mínimo de la función: ', f_rastrigin2d(best_ind[0], best_ind[1]))
print('Variables que producen el valor mínimo: ', best_ind)
Valor mínimo de la función: 0.0 Variables que producen el valor mínimo: (-1.6130253689405727e-09, -1.3999432868862866e-09)
# Definir los parámetros de la evolución diferencial
F = 0.5 # factor de amplificación
CR = 0.9 # tasa de recombinación
NP = 20 # tamaño de la población
ngen = 100 # número de generaciones
# Definir el límite inferior y superior para los valores de x, y, y z
lim_inf = -10
lim_sup = 10
# Generar la población inicial de forma aleatoria
pop = []
for i in range(NP):
x = random.uniform(lim_inf, lim_sup)
y = random.uniform(lim_inf, lim_sup)
z = random.uniform(lim_inf, lim_sup)
pop.append((x, y, z))
# Realizar la evolución diferencial
for gen in range(ngen):
for i in range(NP):
# Seleccionar 3 individuos de forma aleatoria
a, b, c = random.sample(pop, 3)
# Seleccionar un número aleatorio para cada variable
R1, R2, R3 = random.random(), random.random(), random.random()
# Realizar la mutación diferencial
x_new = np.clip(a[0] + F*(b[0] - c[0]), lim_inf, lim_sup)
y_new = np.clip(a[1] + F*(b[1] - c[1]), lim_inf, lim_sup)
z_new = np.clip(a[2] + F*(b[2] - c[2]), lim_inf, lim_sup)
# Realizar la recombinación
if random.random() < CR:
x = x_new
y = y_new
z = z_new
else:
x = a[0]
y = a[1]
z = a[2]
# Evaluar la función en el nuevo individuo
fitness_new = f_rastrigin3d(x, y, z)
# Reemplazar el individuo anterior si el nuevo individuo es mejor
if fitness_new < f_rastrigin3d(pop[i][0], pop[i][1], pop[i][2]):
pop[i] = (x, y, z)
# Imprimir el mejor individuo encontrado
best_ind = min(pop, key=lambda x: f_rastrigin3d(x[0], x[1], x[2]))
print('Valor mínimo de la función: ', f_rastrigin3d(best_ind[0], best_ind[1], best_ind[2]))
print('Variables que producen el valor mínimo: ', best_ind)
Valor mínimo de la función: -5.2820006782696485 Variables que producen el valor mínimo: (0.04069343900956812, 0.09806400831968248, 1.0842032194549434) -5.2820006782696485 Variables que producen el valor mínimo: (0.04069343900956812, 0.09806400831968248, 1.0842032194549434)
def f_griewank2d(x, y):
return 1 + (x**2 + y**2)/4000 - np.cos(x)*np.cos(y/np.sqrt(2))
def f_griewank3d(x, y, z):
return 1 + (x**2 + y**2 + z**2)/4000 - np.cos(x)*np.cos(y/np.sqrt(2))*np.cos(z/np.sqrt(3))
n = 50
x1 = np.linspace(-600, 600, n)
x2 = np.linspace(-600, 600, n)
X1, X2 = np.meshgrid(x1, x2)
X = np.vstack((X1.flatten(), X2.flatten())).T
z = f_griewank2d(X[:, 0], X[:, 1])
Z = z.reshape(n, n)
# Crear un gráfico de contorno con valores numéricos
levels = np.arange(0, 201, 10)
cp = plt.contourf(x1, x2, Z, levels=levels, cmap='plasma')
plt.colorbar(cp)
# Agregar líneas de contorno
plt.contour(x1, x2, Z, levels=levels, colors='k', linewidths=0.5)
# Título y etiquetas de los ejes
plt.title("Función de Griewank")
plt.xlabel("x1")
plt.ylabel("x2")
# Mostrar el gráfico
plt.show()
# Cálculo de la función
z = f_griewank2d(X1, X2)
# Creación de la figura
fig = go.Figure(data=[go.Surface(x=X1, y=X2, z=z)])
# Personalización de la figura
fig.update_layout(
title='Función de Griewank 2D',
scene = dict(
xaxis_title='X1',
yaxis_title='X2',
zaxis_title='griewank(X1,X2)',
aspectmode='cube',
camera=dict(
eye=dict(x=1.75, y=-1.75, z=1.25)
)
)
)
# Mostrar la figura
fig.show()
def griewank(x):
# Esta función calcula el valor de la función Griewank en un punto dado x
n = x.shape[0]
sum = np.sum(x**2)
prod = np.prod(np.cos(x / np.sqrt(np.arange(1, n+1))))
return 1 + sum / 4000 - prod
def griewank_gradient(x):
# Esta función calcula el gradiente de la función Griewank en un punto dado x
n = x.shape[0]
prod = np.prod(np.cos(x / np.sqrt(np.arange(1, n+1))))
grad = np.zeros(n)
for i in range(n):
grad[i] = x[i] / 2000 + (1 / np.sqrt(i+1)) * prod * np.sin(x[i] / np.sqrt(i+1))
return grad
def griewank_hessian(x):
# Esta función calcula la matriz hessiana de la función Griewank en un punto dado x
n = x.shape[0]
prod = np.prod(np.cos(x / np.sqrt(np.arange(1, n+1))))
H = np.zeros((n, n))
for i in range(n):
for j in range(i, n):
if i == j:
H[i][j] = 1 / 2000 + (1 / (i+1)) * prod * (1 - (np.sin(x[i] / np.sqrt(i+1)))**2)
else:
pij = np.prod(np.cos(x / np.sqrt(np.arange(1, n+1)))[:i]) * np.prod(np.cos(x / np.sqrt(np.arange(1, n+1)))[i+1:j])
H[i][j] = (-1 / ((i-j)**2)) * (np.sin(x[i] / np.sqrt(i+1))) * (np.sin(x[j] / np.sqrt(j+1))) * prod * pij
H[j][i] = H[i][j]
return H
def newton_raphson_optimization(x, max_iter=100, eps=1e-6):
# Esta función aplica el método de Newton-Raphson para encontrar el mínimo de la función Griewank
# a partir de un punto inicial x
path = [x]
for i in range(max_iter):
grad = griewank_gradient(x)
H = griewank_hessian(x)
delta_x = np.linalg.solve(H, -grad)
x = x + delta_x
path.append(x)
if np.linalg.norm(delta_x) < eps:
break
return x, path
griewank_scale = 600 / 3 # Scale factor for desired range of [-600, 600] (3 standard deviations)
griewank_init_point_2D = griewank_scale * np.random.randn(2)
griewank_init_point_2D = np.clip(griewank_init_point_2D, -600, 600) # Clip values outside of the desired range
griewank_min_2D, griewank_path_2D = newton_raphson_optimization(griewank_init_point_2D, max_iter= 10000, )
#Impresión de resultados
print(f"Init value: {griewank_init_point_2D}")
print(f"Minimum: {griewank_min_2D}")
print(f"Function value at minimum: {griewank(griewank_min_2D)}")
print("Optimization path: ")
for point in griewank_path_2D:
print(f"Punto evaluado:{point}--", f"Valor de la función en el punto: {griewank(point)}")
Init value: [-268.45690276 -328.4636416 ] Minimum: [-219.79895285 -239.66954303] Function value at minimum: 26.459604466919707 Optimization path: Punto evaluado:[-268.45690276 -328.4636416 ]-- Valor de la función en el punto: 46.1345318918537 Punto evaluado:[-236.8090539 -316.80549964]-- Valor de la función en el punto: 39.89838123232788 Punto evaluado:[-236.92365268 -316.33887784]-- Valor de la función en el punto: 39.838259647444545 Punto evaluado:[-237.09813642 -315.64651031]-- Valor de la función en el punto: 39.871017739133634 Punto evaluado:[-266.2349533 -320.68007028]-- Valor de la función en el punto: 45.019080016296016 Punto evaluado:[-258.46719463 -312.11084895]-- Valor de la función en el punto: 41.590905278617775 Punto evaluado:[-262.4683787 -317.27379128]-- Valor de la función en el punto: 43.427757027481704 Punto evaluado:[-257.43322075 -312.71416805]-- Valor de la función en el punto: 41.66878110225802 Punto evaluado:[-262.42626699 -282.19693939]-- Valor de la función en el punto: 38.120281633940586 Punto evaluado:[-290.61622545 -308.67320533]-- Valor de la función en el punto: 45.932803518253735 Punto evaluado:[-209.4212731 -234.30413624]-- Valor de la función en el punto: 25.36090527066114 Punto evaluado:[-211.75824503 -237.01649732]-- Valor de la función en el punto: 26.11865563006501 Punto evaluado:[-212.02304479 -236.8364216 ]-- Valor de la función en el punto: 26.241675864745016 Punto evaluado:[-219.94341084 -242.45003907]-- Valor de la función en el punto: 28.008635536129148 Punto evaluado:[-220.23055931 -236.56649934]-- Valor de la función en el punto: 27.795748849615975 Punto evaluado:[-219.213345 -240.05651962]-- Valor de la función en el punto: 26.65816286357996 Punto evaluado:[-220.12552779 -239.47943438]-- Valor de la función en el punto: 26.520355625595574 Punto evaluado:[-219.76322165 -239.72024788]-- Valor de la función en el punto: 26.460817235236306 Punto evaluado:[-219.79852135 -239.6708955 ]-- Valor de la función en el punto: 26.459602829915493 Punto evaluado:[-219.79894308 -239.66956946]-- Valor de la función en el punto: 26.45960442541122 Punto evaluado:[-219.79895265 -239.66954355]-- Valor de la función en el punto: 26.459604466104288 Punto evaluado:[-219.79895285 -239.66954303]-- Valor de la función en el punto: 26.459604466919707
# Generate a grid of x, y values
gwk_x_min, gwk_x_max, gwk_x_maxvalues = -600,600, 600
gwk_y_min, gwk_y_max, gwk_y_maxvalues = -600, 600, 600
x_griewank = np.linspace(gwk_x_min, gwk_x_max, gwk_x_maxvalues)
y_griewank = np.linspace(gwk_y_min, gwk_y_max, gwk_y_maxvalues)
X_griewank, Y_griewank = np.meshgrid(x_griewank, y_griewank)
# Compute the Rastrigin function for each (x, y) pair
Z_griewank = np.zeros_like(X_griewank)
for i in range(X_griewank.shape[0]):
for j in range(X_griewank.shape[1]):
Z_griewank[i, j] = griewank(np.array([X_griewank[i, j], Y_griewank[i, j]]))
griewank_minimum= griewank_min_2D.reshape(-1, 1)
#Plot Griewank surface and minimum point founded.
griewank_surface_fig = plt.figure(figsize=(10,8))
griewank_surface_axes = plt.axes(projection='3d', elev=50, azim=-50)
griewank_surface_axes.plot_surface(X_griewank, Y_griewank, Z_griewank, norm=LogNorm(), rstride=1, cstride=1,
edgecolor='none', alpha=1, cmap=plt.cm.jet)
griewank_surface_axes.plot(*griewank_minimum, griewank(griewank_minimum), 'b*', markersize=10)
griewank_surface_axes.set_xlabel('$x$')
griewank_surface_axes.set_ylabel('$y$')
griewank_surface_axes.set_zlabel('$z$')
griewank_surface_axes.set_xlim((gwk_x_min, gwk_x_max))
griewank_surface_axes.set_ylim((gwk_y_min, gwk_y_max))
plt.show()
# Create a 2D contour plot of the Griewank function and plot the optimization path
griewank_contour_fig, griewank_contour_axes = plt.subplots(figsize=(10,8))
griewank_contour = griewank_contour_axes.contour(X_griewank, Y_griewank, Z_griewank, levels=20, cmap='coolwarm')
griewank_contour_axes.set_xlabel('x')
griewank_contour_axes.set_ylabel('y')
for point in griewank_path_2D:
griewank_contour_axes.plot(*point, 'o-',color='blue',alpha=1)
# Add a marker for the global minimum at (0, 0)
griewank_contour_axes.scatter(*griewank_minimum, marker='*', s=500, color='red')
plt.colorbar(griewank_contour)
plt.show()
# Define the function to update the animation
def update(frame):
x = griewank_path_2D[frame]
line.set_data(x[0], x[1])
if frame > 0:
line_trace.set_data([griewank_path_2D[i][0] for i in range(frame+1)], [griewank_path_2D[i][1] for i in range(frame+1)])
return line, line_trace
# Set the starting point for the optimization
# x_start = np.array([-4, 4])
# Run the optimization
# x_opt, path = newton_raphson_optimization(x_start)
# Plot the contour and the starting point
griewank_contour_fig, griewank_contour_axes = plt.subplots(figsize=(10,8))
griewank_contour_axes.set_title('Optimization Path Animation')
griewank_contour_axes.set_xlabel('x')
griewank_contour_axes.set_ylabel('y')
griewank_contour = griewank_contour_axes.contour(X_griewank, Y_griewank, Z_griewank, levels=20, cmap='coolwarm')
griewank_contour_axes.plot(griewank_init_point_2D[0], griewank_init_point_2D[1], 'o', color='blue')
# Set up the animation
line, = griewank_contour_axes.plot([], [], 'o-', color='red', lw=2)
line_trace, = griewank_contour_axes.plot([], [], '-', color='green', lw=4)
ani = animation.FuncAnimation(griewank_contour_fig, update, frames=len(griewank_path_2D), interval=200, blit=True)
HTML(ani.to_jshtml())
griewank_scale = 600 / 3 # Scale factor for desired range of [-600, 600] (3 standard deviations)
griewank_init_point_3D = griewank_scale * np.random.randn(3)
griewank_init_point_3D = np.clip(griewank_init_point_3D, -600, 600) # Clip values outside of the desired range
griewank_min_3D, griewank_path_3D = newton_raphson_optimization(griewank_init_point_3D, max_iter= 10000, )
#Impresión de resultados
print(f"Init value: {griewank_init_point_3D}")
print(f"Minimum: {griewank_min_3D}")
print(f"Function value at minimum: {griewank(griewank_min_3D)}")
print("Optimization path: ")
for point in griewank_path_3D:
print(f"Punto evaluado:{point}--", f"Valor de la función en el punto: {griewank(point)}")
Init value: [ 21.2161394 465.38624143 20.40039014] Minimum: [ 25.11933844 461.55675686 21.7308067 ] Function value at minimum: 53.59734061619003 Optimization path: Punto evaluado:[ 21.2161394 465.38624143 20.40039014]-- Valor de la función en el punto: 55.007908731012144 Punto evaluado:[ 21.15121841 466.11706225 24.21638111]-- Valor de la función en el punto: 55.47495863734799 Punto evaluado:[ 20.70286608 463.18230165 17.23096299]-- Valor de la función en el punto: 54.646545826382294 Punto evaluado:[ 24.06364336 465.14289142 17.82271096]-- Valor de la función en el punto: 55.135178039985966 Punto evaluado:[ 20.76992987 467.28391989 18.25132361]-- Valor de la función en el punto: 55.90864926221509 Punto evaluado:[ 21.51127643 468.4935713 26.83596618]-- Valor de la función en el punto: 56.3079803349592 Punto evaluado:[ 25.10407543 461.50735351 22.57318044]-- Valor de la función en el punto: 53.70658272036371 Punto evaluado:[ 25.12726044 462.02720743 21.21429798]-- Valor de la función en el punto: 53.68813644261122 Punto evaluado:[ 25.1196935 461.58142364 21.76634755]-- Valor de la función en el punto: 53.59729981576747 Punto evaluado:[ 25.11932819 461.55895552 21.73104426]-- Valor de la función en el punto: 53.597307869603654 Punto evaluado:[ 25.11933676 461.55693182 21.73078482]-- Valor de la función en el punto: 53.59733792691605 Punto evaluado:[ 25.1193383 461.5567711 21.73080493]-- Valor de la función en el punto: 53.597340396759385 Punto evaluado:[ 25.11933843 461.55675802 21.73080656]-- Valor de la función en el punto: 53.5973405984328 Punto evaluado:[ 25.11933844 461.55675695 21.73080669]-- Valor de la función en el punto: 53.59734061485335 Punto evaluado:[ 25.11933844 461.55675686 21.7308067 ]-- Valor de la función en el punto: 53.59734061619003
# Definir los límites de los parámetros
limite_inferior = -5.0
limite_superior = 5.0
# Definir los parámetros del algoritmo genético
num_generaciones = 50
tamano_poblacion = 100
probabilidad_mutacion = 0.2
probabilidad_cruce = 0.5
# Crear los tipos de objeto para el problema
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individuo", list, fitness=creator.FitnessMin)
# Crear la caja de herramientas del algoritmo genético
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, limite_inferior, limite_superior)
toolbox.register("individuo", tools.initCycle, creator.Individuo,
(toolbox.attr_float, toolbox.attr_float), n=1)
toolbox.register("poblacion", tools.initRepeat, list, toolbox.individuo)
# Definir la función de evaluación
def evaluar(individuo):
return f_griewank2d(individuo[0], individuo[1]),
# Definir los operadores genéticos
toolbox.register("evaluate", evaluar)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)
# Crear la población inicial
poblacion = toolbox.poblacion(n=tamano_poblacion)
# Ejecutar el algoritmo genético
for generacion in range(num_generaciones):
# Seleccionar a los padres
padres = toolbox.select(poblacion, len(poblacion))
# Clonar los padres para aplicar operadores genéticos
hijos = list(map(toolbox.clone, padres))
# Aplicar cruce y mutación a los hijos
for hijo in hijos:
if random.random() < probabilidad_cruce:
toolbox.mate(hijo, random.choice(hijos))
del hijo.fitness.values
if random.random() < probabilidad_mutacion:
toolbox.mutate(hijo)
del hijo.fitness.values
# Evaluar la aptitud de los nuevos individuos
individuos_nuevos = [individuo for individuo in hijos if not individuo.fitness.valid]
aptitudes_nuevas = toolbox.map(toolbox.evaluate, individuos_nuevos)
for individuo, aptitud in zip(individuos_nuevos, aptitudes_nuevas):
individuo.fitness.values = aptitud
# Reemplazar la población anterior por la nueva generación
poblacion = hijos
# Seleccionar al mejor individuo
mejor_individuo = tools.selBest(poblacion, 1)[0]
print("Mejor solución encontrada:", mejor_individuo)
print("Valor de la función objetivo:", mejor_individuo.fitness.values[0])
Mejor solución encontrada: [3.3199442668675464e-06, -2.4718480383282513e-06] Valor de la función objetivo: 7.042810779012143e-12
c:\Users\Hp\anaconda3\envs\tf_gpu\lib\site-packages\deap\creator.py:138: RuntimeWarning: A class named 'FitnessMin' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it. c:\Users\Hp\anaconda3\envs\tf_gpu\lib\site-packages\deap\creator.py:138: RuntimeWarning: A class named 'Individuo' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it.
# Definir los límites de los parámetros
limite_inferior = -5.0
limite_superior = 5.0
# Definir los parámetros del algoritmo genético
num_generaciones = 50
tamano_poblacion = 100
probabilidad_mutacion = 0.2
probabilidad_cruce = 0.5
# Crear los tipos de objeto para el problema
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individuo", list, fitness=creator.FitnessMin)
# Crear la caja de herramientas del algoritmo genético
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, limite_inferior, limite_superior)
toolbox.register("individuo", tools.initRepeat, creator.Individuo,
toolbox.attr_float, n=3)
toolbox.register("poblacion", tools.initRepeat, list, toolbox.individuo)
# Definir la función de evaluación
def evaluar(individuo):
return f_griewank3d(individuo[0], individuo[1], individuo[2]),
# Definir los operadores genéticos
toolbox.register("evaluate", evaluar)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)
# Crear la población inicial
poblacion = toolbox.poblacion(n=tamano_poblacion)
# Ejecutar el algoritmo genético
for generacion in range(num_generaciones):
# Seleccionar a los padres
padres = toolbox.select(poblacion, len(poblacion))
# Clonar los padres para aplicar operadores genéticos
hijos = list(map(toolbox.clone, padres))
# Aplicar cruce y mutación a los hijos
for hijo in hijos:
if random.random() < probabilidad_cruce:
toolbox.mate(hijo, random.choice(hijos))
del hijo.fitness.values
if random.random() < probabilidad_mutacion:
toolbox.mutate(hijo)
del hijo.fitness.values
# Evaluar la aptitud de los nuevos individuos
individuos_nuevos = [individuo for individuo in hijos if not individuo.fitness.valid]
aptitudes_nuevas = toolbox.map(toolbox.evaluate, individuos_nuevos)
for individuo, aptitud in zip(individuos_nuevos, aptitudes_nuevas):
individuo.fitness.values = aptitud
# Reemplazar la población anterior por la nueva generación
poblacion = hijos
# Seleccionar al mejor individuo
mejor_individuo = tools.selBest(poblacion, 1)[0]
print("Mejor solución encontrada:", mejor_individuo)
print("Valor de la función objetivo:", mejor_individuo.fitness.values[0])
Mejor solución encontrada: [3.307304896771483e-07, 1.4990202736401472e-06, -6.790513197076734e-07] Valor de la función objetivo: 6.940004126931854e-13
#Definimos la funcion Griewank para los calculos
def f(x, y):
return (x ** 2/(4000) + y ** 2/(4000)) - ((np.cos(x/np.sqrt(1))+1) * (np.cos(y/np.sqrt(2))+1))
#Calculamos los datos X Y de la funcion
x, y = np.array(np.meshgrid(np.linspace(-600, 600, 100),
np.linspace(-600, 600, 100)))
#Evalamos los datos X Y
z = f(x, y)
#Encontramos el minimo Global de la funcion
x_min = x.ravel()[z.argmin()]
y_min = y.ravel()[z.argmin()]
# Se inicializan los Hiperparametros del algoritmo PSO
c1 = c2 = 0.1 # Factores de Acelearacion
w = 0.8 # Coeficiente de Inercia
# Creamos las particulas
n_particles = 35
np.random.seed(100)
X = np.random.rand(2, n_particles) * 500 # Posiscion Inicial de las particulas
V = np.random.randn(2, n_particles) * 200 # velocidad de las particulas
# Inicializacion de datos
pbest = X # Mejor Posicion de la particula
pbest_obj = f(X[0], X[1])
gbest = pbest[:, pbest_obj.argmin()] # Mejor Posicion del Enjambre
gbest_obj = pbest_obj.min()
#funciond de actualizacion de los puntos del enjambre
def update():
global V, X, pbest, pbest_obj, gbest, gbest_obj
# actualizacion de parametros
r1, r2 = np.random.rand(2)
# Actualizacion de velocidades
V = w * V + c1*r1*(pbest - X) + c2*r2*(gbest.reshape(-1, 1)-X)
X = X + V # Actualizacion de Posicion de particulas
obj = f(X[0], X[1]) # Nueva Posicion
# Se evalua si la posicion nueva 'obj' es mejor a la anterior
pbest[:, (pbest_obj >= obj)] = X[:, (pbest_obj >= obj)]
pbest_obj = np.array([pbest_obj, obj]).min(axis=0)
gbest = pbest[:, pbest_obj.argmin()]
gbest_obj = pbest_obj.min()
# Visualizacion del grafico de Grienwank Y los puntos del enjambre
fig, ax = plt.subplots(figsize=(16, 8))
fig.set_tight_layout(True)
img = ax.imshow(z, extent=[-600, 600, -600, 600],
origin='lower', cmap='plasma', alpha=0.5)
fig.colorbar(img, ax=ax)
ax.plot([x_min], [y_min], marker='x', markersize=10,
color="white") # Minimo Global
contours = ax.contour(x, y, z, 10, colors='black', alpha=0.4)
ax.clabel(contours, inline=True, fontsize=12, fmt="%.0f")
# mejor posision del punto
pbest_plot = ax.scatter(
pbest[0], pbest[1], marker='o', color='black', alpha=0.5)
p_plot = ax.scatter(X[0], X[1], marker='o', color='green', alpha=0.5) # Puntos
p_arrow = ax.quiver(X[0], X[1], V[0], V[1], color='green', width=0.005,
angles='xy', scale_units='xy', scale=1) # flecha de los puntos
gbest_plot = plt.scatter([gbest[0]], [gbest[1]], marker='*',
s=100, color='red', alpha=0.4) # Mejor posicion del enjambre
ax.set_xlim([-600, 600])
ax.set_ylim([-600, 600])
def animate(i):
# Animacion de los puntos en moviento
title = 'Iteration {:02d}'.format(i)
# Update params
update()
# Set picture
ax.set_title(title)
pbest_plot.set_offsets(pbest.T)
p_plot.set_offsets(X.T)
p_arrow.set_offsets(X.T)
p_arrow.set_UVC(V[0], V[1])
gbest_plot.set_offsets(gbest.reshape(1, -1))
return ax, pbest_plot, p_plot, p_arrow, gbest_plot
anim = FuncAnimation(fig, animate, frames=list(
range(1, 50)), interval=500, blit=False, repeat=True)
# Se guarda la animacion en formato gif
anim.save("PSOGRI.gif", dpi=120, writer="imagemagick")
# Informacion genreal del algoritmo luego de optimizar
print("PSO found best solution at f({})={}".format(gbest, gbest_obj))
print("Global optimal at f({})={}".format([x_min,y_min], f(x_min,y_min)))
MovieWriter imagemagick unavailable; using Pillow instead.
PSO found best solution at f([-6.33640917 0.14168459])=-3.977103581529073 Global optimal at f([-6.060606060606119, 18.18181818181813])=-3.776287411893872
from IPython.display import Image, display
Image(filename='PSOGRI.gif')
<IPython.core.display.Image object>
lb = [-10, -10, -10] # límite inferior de las variables
ub = [10, 10, 10] # límite superior de las variables
x_opt, f_opt = pso(lambda x: f_griewank3d(x[0], x[1], x[2]), lb, ub)
print("El punto óptimo es:", x_opt)
Stopping search: Swarm best objective change less than 1e-08 El punto óptimo es: [-6.27996580e+00 5.09301382e-05 7.34195399e-05]
# Definir los parámetros de la evolución diferencial
F = 0.5 # factor de amplificación
CR = 0.9 # tasa de recombinación
NP = 20 # tamaño de la población
ngen = 100 # número de generaciones
# Definir el límite inferior y superior para los valores de x e y
lim_inf = -10
lim_sup = 10
# Generar la población inicial de forma aleatoria
pop = []
for i in range(NP):
x = random.uniform(lim_inf, lim_sup)
y = random.uniform(lim_inf, lim_sup)
pop.append((x, y))
# Realizar la evolución diferencial
for gen in range(ngen):
for i in range(NP):
# Seleccionar 3 individuos de forma aleatoria
a, b, c = random.sample(pop, 3)
# Seleccionar un número aleatorio para cada variable
R1, R2 = random.random(), random.random()
# Realizar la mutación diferencial
x_new = np.clip(a[0] + F*(b[0] - c[0]), lim_inf, lim_sup)
y_new = np.clip(a[1] + F*(b[1] - c[1]), lim_inf, lim_sup)
# Realizar la recombinación
if random.random() < CR:
x = x_new
y = y_new
else:
x = a[0]
y = a[1]
# Evaluar la función en el nuevo individuo
fitness_new = f_griewank2d(x, y)
# Reemplazar el individuo anterior si el nuevo individuo es mejor
if fitness_new < f_griewank2d(pop[i][0], pop[i][1]):
pop[i] = (x, y)
# Imprimir el mejor individuo encontrado
best_ind = min(pop, key=lambda x: f_griewank2d(x[0], x[1]))
print('Valor mínimo de la función: ', f_griewank2d(best_ind[0], best_ind[1]))
print('Variables que producen el valor mínimo: ', best_ind)
Valor mínimo de la función: 0.007396040334114895 Variables que producen el valor mínimo: (3.140022634806684, 4.438444483718577)
# Definir los parámetros de la evolución diferencial
F = 0.5 # factor de amplificación
CR = 0.9 # tasa de recombinación
NP = 20 # tamaño de la población
ngen = 100 # número de generaciones
# Definir el límite inferior y superior para los valores de x, y, y z
lim_inf = -10
lim_sup = 10
# Generar la población inicial de forma aleatoria
pop = []
for i in range(NP):
x = random.uniform(lim_inf, lim_sup)
y = random.uniform(lim_inf, lim_sup)
z = random.uniform(lim_inf, lim_sup)
pop.append((x, y, z))
# Realizar la evolución diferencial
for gen in range(ngen):
for i in range(NP):
# Seleccionar 3 individuos de forma aleatoria
a, b, c = random.sample(pop, 3)
# Seleccionar un número aleatorio para cada variable
R1, R2, R3 = random.random(), random.random(), random.random()
# Realizar la mutación diferencial
x_new = np.clip(a[0] + F*(b[0] - c[0]), lim_inf, lim_sup)
y_new = np.clip(a[1] + F*(b[1] - c[1]), lim_inf, lim_sup)
z_new = np.clip(a[2] + F*(b[2] - c[2]), lim_inf, lim_sup)
# Realizar la recombinación
if random.random() < CR:
x = x_new
y = y_new
z = z_new
else:
x = a[0]
y = a[1]
z = a[2]
# Evaluar la función en el nuevo individuo
fitness_new = f_griewank3d(x, y, z)
# Reemplazar el individuo anterior si el nuevo individuo es mejor
if fitness_new < f_griewank3d(pop[i][0], pop[i][1], pop[i][2]):
pop[i] = (x, y, z)
# Imprimir el mejor individuo encontrado
best_ind = min(pop, key=lambda x: f_griewank3d(x[0], x[1], x[2]))
print('Valor mínimo de la función: ', f_griewank3d(best_ind[0], best_ind[1], best_ind[2]))
print('Variables que producen el valor mínimo: ', best_ind)
Valor mínimo de la función: 0.056348699247850464 Variables que producen el valor mínimo: (3.062545509753079, 0.38435240630277434, -5.639937371581347)
Parte 2.1: Problema Viajero Colombian Version.¶
#Importamos las librerias Necesarias para el desarollo Utlizando Algoritmos Geneticos
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import Point
import pyproj
from pygad import gann
from matplotlib import animation, rc
import pandas as pd
import random
import itertools
import imageio
import os
import numpy as np
import matplotlib.pyplot as plt
import time
import warnings
warnings.filterwarnings("ignore")
class AntColonyOptimizer:
def __init__(self, ants, evaporation_rate, intensification, alpha=1.0, beta=0.0, beta_evaporation_rate=0,
choose_best=.1):
"""
Ant colony optimizer. Traverses a graph and finds either the max or min distance between nodes.
:param ants: number of ants to traverse the graph
:param evaporation_rate: rate at which pheromone evaporates
:param intensification: constant added to the best path
:param alpha: weighting of pheromone
:param beta: weighting of heuristic (1/distance)
:param beta_evaporation_rate: rate at which beta decays (optional)
:param choose_best: probability to choose the best route
"""
# Parameters
self.ants = ants
self.evaporation_rate = evaporation_rate
self.pheromone_intensification = intensification
self.heuristic_alpha = alpha
self.heuristic_beta = beta
self.beta_evaporation_rate = beta_evaporation_rate
self.choose_best = choose_best
# Internal representations
self.pheromone_matrix = None
self.heuristic_matrix = None
self.probability_matrix = None
self.map = None
self.set_of_available_nodes = None
# Internal stats
self.best_series = []
self.best = None
self.fitted = False
self.best_path = None
self.fit_time = None
# Plotting values
self.stopped_early = False
def __str__(self):
string = "Ant Colony Optimizer"
string += "\n--------------------"
string += "\nDesigned to optimize either the minimum or maximum distance between nodes in a square matrix that behaves like a distance matrix."
string += "\n--------------------"
string += "\nNumber of ants:\t\t\t\t{}".format(self.ants)
string += "\nEvaporation rate:\t\t\t{}".format(self.evaporation_rate)
string += "\nIntensification factor:\t\t{}".format(
self.pheromone_intensification)
string += "\nAlpha Heuristic:\t\t\t{}".format(self.heuristic_alpha)
string += "\nBeta Heuristic:\t\t\t\t{}".format(self.heuristic_beta)
string += "\nBeta Evaporation Rate:\t\t{}".format(
self.beta_evaporation_rate)
string += "\nChoose Best Percentage:\t\t{}".format(self.choose_best)
string += "\n--------------------"
string += "\nUSAGE:"
string += "\nNumber of ants influences how many paths are explored each iteration."
string += "\nThe alpha and beta heuristics affect how much influence the pheromones or the distance heuristic weigh an ants' decisions."
string += "\nBeta evaporation reduces the influence of the heuristic over time."
string += "\nChoose best is a percentage of how often an ant will choose the best route over probabilistically choosing a route based on pheromones."
string += "\n--------------------"
if self.fitted:
string += "\n\nThis optimizer has been fitted."
else:
string += "\n\nThis optimizer has NOT been fitted."
return string
def _initialize(self):
"""
Initializes the model by creating the various matrices and generating the list of available nodes
"""
assert self.map.shape[0] == self.map.shape[1], "Map is not a distance matrix!"
num_nodes = self.map.shape[0]
self.pheromone_matrix = np.ones((num_nodes, num_nodes))
# Remove the diagonal since there is no pheromone from node i to itself
self.pheromone_matrix[np.eye(num_nodes) == 1] = 0
self.heuristic_matrix = 1 / self.map
self.probability_matrix = (self.pheromone_matrix ** self.heuristic_alpha) * (
self.heuristic_matrix ** self.heuristic_beta) # element by element multiplcation
self.set_of_available_nodes = list(range(num_nodes))
def _reinstate_nodes(self):
"""
Resets available nodes to all nodes for the next iteration
"""
self.set_of_available_nodes = list(range(self.map.shape[0]))
def _update_probabilities(self):
"""
After evaporation and intensification, the probability matrix needs to be updated. This function
does that.
"""
self.probability_matrix = (self.pheromone_matrix ** self.heuristic_alpha) * (
self.heuristic_matrix ** self.heuristic_beta)
def _choose_next_node(self, from_node):
"""
Chooses the next node based on probabilities. If p < p_choose_best, then the best path is chosen, otherwise
it is selected from a probability distribution weighted by the pheromone.
:param from_node: the node the ant is coming from
:return: index of the node the ant is going to
"""
numerator = self.probability_matrix[from_node,
self.set_of_available_nodes]
if np.random.random() < self.choose_best:
next_node = np.argmax(numerator)
else:
denominator = np.sum(numerator)
probabilities = numerator / denominator
next_node = np.random.choice(
range(len(probabilities)), p=probabilities)
return next_node
def _remove_node(self, node):
self.set_of_available_nodes.remove(node)
def _evaluate(self, paths, mode):
"""
Evaluates the solutions of the ants by adding up the distances between nodes.
:param paths: solutions from the ants
:param mode: max or min
:return: x and y coordinates of the best path as a tuple, the best path, and the best score
"""
scores = np.zeros(len(paths))
coordinates_i = []
coordinates_j = []
for index, path in enumerate(paths):
score = 0
coords_i = []
coords_j = []
for i in range(len(path) - 1):
coords_i.append(path[i])
coords_j.append(path[i + 1])
score += self.map[path[i], path[i + 1]]
scores[index] = score
coordinates_i.append(coords_i)
coordinates_j.append(coords_j)
if mode == 'min':
best = np.argmin(scores)
elif mode == 'max':
best = np.argmax(scores)
return (coordinates_i[best], coordinates_j[best]), paths[best], scores[best]
def _evaporation(self):
"""
Evaporate some pheromone as the inverse of the evaporation rate. Also evaporates beta if desired.
"""
self.pheromone_matrix *= (1 - self.evaporation_rate)
self.heuristic_beta *= (1 - self.beta_evaporation_rate)
def _intensify(self, best_coords):
"""
Increases the pheromone by some scalar for the best route.
:param best_coords: x and y (i and j) coordinates of the best route
"""
i = best_coords[0]
j = best_coords[1]
self.pheromone_matrix[i, j] += self.pheromone_intensification
def fit(self, map_matrix, iterations=100, mode='min', early_stopping_count=20, verbose=True):
"""
Fits the ACO to a specific map. This was designed with the Traveling Salesman problem in mind.
:param map_matrix: Distance matrix or some other matrix with similar properties
:param iterations: number of iterations
:param mode: whether to get the minimum path or maximum path
:param early_stopping_count: how many iterations of the same score to make the algorithm stop early
:return: the best score
"""
if verbose:
print("Beginning ACO Optimization with {} iterations...".format(iterations))
self.map = map_matrix
start = time.time()
self._initialize()
num_equal = 0
for i in range(iterations):
start_iter = time.time()
paths = []
path = []
for ant in range(self.ants):
current_node = self.set_of_available_nodes[np.random.randint(
0, len(self.set_of_available_nodes))]
start_node = current_node
while True:
path.append(current_node)
self._remove_node(current_node)
if len(self.set_of_available_nodes) != 0:
current_node_index = self._choose_next_node(
current_node)
current_node = self.set_of_available_nodes[current_node_index]
else:
break
path.append(start_node) # go back to start
self._reinstate_nodes()
paths.append(path)
path = []
best_path_coords, best_path, best_score = self._evaluate(
paths, mode)
if i == 0:
best_score_so_far = best_score
else:
if mode == 'min':
if best_score < best_score_so_far:
best_score_so_far = best_score
self.best_path = best_path
elif mode == 'max':
if best_score > best_score_so_far:
best_score_so_far = best_score
self.best_path = best_path
if best_score == best_score_so_far:
num_equal += 1
else:
num_equal = 0
self.best_series.append(best_score)
self._evaporation()
self._intensify(best_path_coords)
self._update_probabilities()
if verbose:
print("Best score at iteration {}: {}; overall: {} ({}s)"
"".format(i, round(best_score, 2), round(best_score_so_far, 2),
round(time.time() - start_iter)))
if best_score == best_score_so_far and num_equal == early_stopping_count:
self.stopped_early = True
print("Stopping early due to {} iterations of the same score.".format(
early_stopping_count))
break
self.fit_time = round(time.time() - start)
self.fitted = True
if mode == 'min':
self.best = self.best_series[np.argmin(self.best_series)]
if verbose:
print(
"ACO fitted. Runtime: {} minutes. Best score: {}".format(self.fit_time // 60, self.best))
return self.best
elif mode == 'max':
self.best = self.best_series[np.argmax(self.best_series)]
if verbose:
print(
"ACO fitted. Runtime: {} minutes. Best score: {}".format(self.fit_time // 60, self.best))
return self.best
else:
raise ValueError("Invalid mode! Choose 'min' or 'max'.")
def plot(self):
"""
Plots the score over time after the model has been fitted.
:return: None if the model isn't fitted yet
"""
if not self.fitted:
print("Ant Colony Optimizer not fitted! There exists nothing to plot.")
return None
else:
fig, ax = plt.subplots(figsize=(20, 15))
ax.plot(self.best_series, label="Best Run")
ax.set_xlabel("Iteration")
ax.set_ylabel("Performance")
ax.text(.8, .6,
'Ants: {}\nEvap Rate: {}\nIntensify: {}\nAlpha: {}\nBeta: {}\nBeta Evap: {}\nChoose Best: {}\n\nFit Time: {}m{}'.format(
self.ants, self.evaporation_rate, self.pheromone_intensification, self.heuristic_alpha,
self.heuristic_beta, self.beta_evaporation_rate, self.choose_best, self.fit_time // 60,
["\nStopped Early!" if self.stopped_early else ""][0]),
bbox={'facecolor': 'gray', 'alpha': 0.8, 'pad': 10}, transform=ax.transAxes)
ax.legend()
plt.title("Ant Colony Optimization Results (best: {})".format(
np.round(self.best, 2)))
plt.show()
Instanciamos las Ciudades y las coordenadas
# Coordenades X Y de las ciudades
ciudades = ["Palmira", "Pasto", "Tuluá", "Bogota", "Pereira", "Armenia", "Manizales", "Valledupar",
"Montería", "Soledad", "Cartagena", "Barranquilla", "Medellín", "Bucaramanga", "Cúcuta"]
idx_ciudades=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14]
#Las coordenadas corresponden al orden en como esta la lista de ciudades
coords=[[3.53944, -76.30361], [1.21361, -77.28111], [4.08466, -76.19536], [4.60971, -74.08175], [4.81333, -75.69611], [4.53389, -75.68111], [5.06889, -75.51738], [10.46314, -73.25322],
[8.74798, -75.88143], [10.91843, -74.76459], [10.39972, -75.51444], [10.96854, -74.78132], [6.25184, -75.56359], [7.12539, -73.1198], [7.89391, -72.50782]]
Y = [coords[i][0] for i in range(len(coords))]
X = [coords[i][1] for i in range(len(coords))]
print(X)
print(Y)
[-76.30361, -77.28111, -76.19536, -74.08175, -75.69611, -75.68111, -75.51738, -73.25322, -75.88143, -74.76459, -75.51444, -74.78132, -75.56359, -73.1198, -72.50782] [3.53944, 1.21361, 4.08466, 4.60971, 4.81333, 4.53389, 5.06889, 10.46314, 8.74798, 10.91843, 10.39972, 10.96854, 6.25184, 7.12539, 7.89391]
Se Realiza una conversion de loas coordenadas Anteriores para que coincidan con el formato del archvo .shp del mapa de Colombia
#TRANSFORMACION DE LAS COORDENAS
# Definir el sistema de proyección de origen (WGS84)
source_crs = 'EPSG:4326'
# Definir el sistema de proyección de destino
target_crs = 'EPSG:21897'
# Crear un objeto de transformación de proyección
transformer = pyproj.Transformer.from_crs(
source_crs, target_crs, always_xy=True)
X_T, Y_T = transformer.transform(X, Y)
print(X_T)
print(Y_T)
# Crear una lista de puntos Shapely a partir de las coordenadas X e Y
geometry = [Point(x, y) for x, y in zip(X_T, Y_T)]
# Crear un GeoDataFrame con las ciudades y la geometría
df = gpd.GeoDataFrame(
{'Ciudad': ciudades, 'geometry': geometry})
#Cambia el formato de POINT a FLOAT puesto que mas adelanto realizaremos calculos con estos
def punto_a_entero(punto):
x, y = punto.x, punto.y
return float(x), float(y)
#Agregamos la columna Punto_float al dataframe
df['punto_float'] = df['geometry'].apply(lambda punto: punto_a_entero(punto))
df
[752605.9915458227, 643275.0465845948, 764787.6710086607, 999529.3630757327, 820426.0323943375, 822021.5505042156, 840318.4461557847, 1090247.024864841, 801474.1002934285, 924887.5034194184, 842633.5437126327, 923071.5071483594, 835526.5527440282, 1105793.363446178, 1173108.296118954] [883431.2205430069, 626167.8831370582, 943734.5434153986, 1001493.9629269837, 1024217.0013714742, 993300.2959815342, 1052441.7478454192, 1648963.4840482422, 1459613.1832835586, 1699295.6026588602, 1642191.7124144642, 1704843.1924954015, 1183313.114715379, 1279803.6634820616, 1365014.5786278923]
| Ciudad | geometry | punto_float | |
|---|---|---|---|
| 0 | Palmira | POINT (752605.992 883431.221) | (752605.9915458227, 883431.2205430069) |
| 1 | Pasto | POINT (643275.047 626167.883) | (643275.0465845948, 626167.8831370582) |
| 2 | Tuluá | POINT (764787.671 943734.543) | (764787.6710086607, 943734.5434153986) |
| 3 | Bogota | POINT (999529.363 1001493.963) | (999529.3630757327, 1001493.9629269837) |
| 4 | Pereira | POINT (820426.032 1024217.001) | (820426.0323943375, 1024217.0013714742) |
| 5 | Armenia | POINT (822021.551 993300.296) | (822021.5505042156, 993300.2959815342) |
| 6 | Manizales | POINT (840318.446 1052441.748) | (840318.4461557847, 1052441.7478454192) |
| 7 | Valledupar | POINT (1090247.025 1648963.484) | (1090247.024864841, 1648963.4840482422) |
| 8 | Montería | POINT (801474.100 1459613.183) | (801474.1002934285, 1459613.1832835586) |
| 9 | Soledad | POINT (924887.503 1699295.603) | (924887.5034194184, 1699295.6026588602) |
| 10 | Cartagena | POINT (842633.544 1642191.712) | (842633.5437126327, 1642191.7124144642) |
| 11 | Barranquilla | POINT (923071.507 1704843.192) | (923071.5071483594, 1704843.1924954015) |
| 12 | Medellín | POINT (835526.553 1183313.115) | (835526.5527440282, 1183313.114715379) |
| 13 | Bucaramanga | POINT (1105793.363 1279803.663) | (1105793.363446178, 1279803.6634820616) |
| 14 | Cúcuta | POINT (1173108.296 1365014.579) | (1173108.296118954, 1365014.5786278923) |
# Cargar el shapefile de Colombia
mapa = gpd.read_file('Colombia Shape\depto.shp')
#print(mapa.crs) #Formato CRS del archivo .shp
fig, ax = plt.subplots(figsize=(14, 14))
mapa.plot(ax=ax, color='lime', edgecolor='black')
df.plot(ax=ax, markersize=50, color='red')
for i, txt in enumerate(df['Ciudad']):
ax.annotate(txt, (df['geometry'][i].x,df['geometry'][i].y), fontsize=12)
Matrices Requeridas:
MATRIZ DE DISTANCIAS ENTRE LAS CIUDADES:MDC
MATRIZ DE # DE PEJAES ENTRE LAS CIUDADES:MNP
MATRIZ DE TIEMPO ESTIMADO ENTRE CIUDADES:MTC
Para realizar estos calculos tambien hay que tener en cuenta la informacion del Automovil que se empleara para este problema el cual es un Chevrolet Camaro SS Este Camaro usa Gasolina como combustible (Inserte Imagen de Camaro)
De este Auto tenemos lo siguiente:
Su rendimeinto en KM por Galon es: 29.15Km/Galon.
El precio de la gasolina a la fecha de realizacion de este trabajo (Marzo 15 del 2023) es de $9.664 Pesos Colombianos por Galon.
EL sueldo por Hora de nuestro chofer a la fecha de realizacion de este trabajo (Marzo 15 del 2023) sera de $1368,39 COP
Nuestro automovil en Colombia pagaria la tarifa de categoria 1. Este costo sera de Aproximadamente $8000 COP
Para esto Usaremos la Formula de Haversine que sirve para calcular la distancia en km entre 2 puntos de la tierra el cual nos dara un estimado de las diatancias que requerimos para la solucion de este problema.
#FUNCION DE HAVERSINE
from math import radians, sin, cos, sqrt, atan2
def haversine(coord1, coord2):
# Radio de la Tierra en km
R = 6371.0
# Convertir coordenadas de grados a radianes
lat1, lon1 = map(radians, coord1)
lat2, lon2 = map(radians, coord2)
# Diferencia de latitud y longitud
dlat = lat2 - lat1
dlon = lon2 - lon1
# Fórmula de Haversine
a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))
# Distancia entre las dos coordenadas
distance = R * c
return distance
# Creamos una matriz vacía de distancias
distances = [[0.0] * len(coords) for i in range(len(coords))]
# Calcular la distancia entre cada par de ciudades
for i in range(len(coords)):
for j in range(i + 1, len(coords)):
# Calcular la distancia entre las coordenadas
dist = haversine(coords[i], coords[j])
# Almacenar la distancia en la matriz y los redondeamos a 2 decimales
distances[i][j] = round(dist,2)
distances[j][i] = round(dist,2)
# Imprimir la matriz de distancias
#for i in range(len(ciudades)):
# print(ciudades[i], distances[i])
#convertimos la lista en una matriz 15X15
MDC = np.array(distances).reshape(15,15)
MDC
array([[ 0. , 280.49, 61.8 , 273.66, 156.85, 130.36, 191.11,
840.18, 581.04, 837.86, 767.78, 842.95, 312.55, 532.17,
640.89],
[ 280.49, 0. , 341.26, 518.46, 437.24, 409.72, 471.28,
1120.66, 851.99, 1114.34, 1039.94, 1119.27, 591.73, 803.06,
911.8 ],
[ 61.8 , 341.26, 0. , 241.51, 98.12, 75.81, 132.76,
779.9 , 519.7 , 776.06, 706.2 , 781.14, 250.93, 479.72,
587.92],
[ 273.66, 518.46, 241.51, 0. , 180.33, 177.48, 167.06,
657.24, 501.23, 705.52, 662.89, 711.25, 245.45, 299.28,
404.5 ],
[ 156.85, 437.24, 98.12, 180.33, 0. , 31.12, 34.63,
683.44, 437.99, 686.56, 621.5 , 691.8 , 160.63, 383.75,
491.38],
[ 130.36, 409.72, 75.81, 177.48, 31.12, 0. , 62.19,
711.52, 469.11, 717.07, 652.51, 722.33, 191.47, 404.1 ,
512.45],
[ 191.11, 471.28, 132.76, 167.06, 34.63, 62.19, 0. ,
649.58, 411.06, 655.7 , 592.76, 660.99, 131.64, 350.08,
457.4 ],
[ 840.18, 1120.66, 779.9 , 657.24, 683.44, 711.52, 649.58,
0. , 345.53, 172.72, 247.38, 176.16, 532.78, 371.43,
297.17],
[ 581.04, 851.99, 519.7 , 501.23, 437.99, 469.11, 411.06,
345.53, 0. , 270.59, 188.02, 274.75, 279.76, 353.62,
383.13],
[ 837.86, 1114.34, 776.06, 705.52, 686.56, 717.07, 655.7 ,
172.72, 270.59, 0. , 100.2 , 5.86, 526.28, 458.8 ,
417.59],
[ 767.78, 1039.94, 706.2 , 662.89, 621.5 , 652.51, 592.76,
247.38, 188.02, 100.2 , 0. , 102.07, 461.25, 449.22,
431.93],
[ 842.95, 1119.27, 781.14, 711.25, 691.8 , 722.33, 660.99,
176.16, 274.75, 5.86, 102.07, 0. , 531.47, 464.64,
423.16],
[ 312.55, 591.73, 250.93, 245.45, 160.63, 191.47, 131.64,
532.78, 279.76, 526.28, 461.25, 531.47, 0. , 286.83,
383.45],
[ 532.17, 803.06, 479.72, 299.28, 383.75, 404.1 , 350.08,
371.43, 353.62, 458.8 , 449.22, 464.64, 286.83, 0. ,
108.88],
[ 640.89, 911.8 , 587.92, 404.5 , 491.38, 512.45, 457.4 ,
297.17, 383.13, 417.59, 431.93, 423.16, 383.45, 108.88,
0. ]])
Ahora Procedemos a las creacion de las matrices MNP y MTC las que corresponden al numero de peajes entre ciudades y el tiempo para viajar de una ciudad a otra
#procedemos a convertir los archivos .xlsx a matrices
# Cargar el archivo de Excel en un DataFrame
Peajes = pd.read_excel('PeajesSimples.xlsx',header=None)
# Convertir el DataFrame en un array de NumPy
MNP = Peajes.to_numpy()
Tiempos = pd.read_excel('TiemposSimples.xlsx', header=None) #en Total de Minutos
# Convertir el DataFrame en un array de NumPy
MTC = Tiempos.to_numpy()
print(MNP.shape)
print(MTC.shape)
(15, 15) (15, 15)
Con estas 3 matrices ya calculadas damos luz para la matriz de costos totales de movilizacion entre ciudades la denotaremos MCTC
#Datos usados para los caclculos
Vh = 1368.39 # Salrio de una hora de sueldo del chofer
Cpc1 = 8000 # Costo Aproximado de los peajes (Esto puede variar mucho o poco entre peajes para la practicidad del ejercicio se fijo).
Pg = 9.664 # PRecio del Galon de Combustible
#Realizamos alguno Ajustes en las matrices
Txh = MTC/60 # Convertimos los minutos en horas
Dxg = MDC/29.15 # Con esta operacion conseguimos el numero de kilometros por galon
#FINALMENTE Calculamos la matriz MCTC
MCTC = Vh*Txh + MNP*8000 + Pg*Dxg
print(MCTC.shape)
(15, 15)
problem = MCTC # matriz de costos totales
Hormigas=[10,20,80,100,120] # nuemro de Hormigas a evaluar
MejorRuta=[] #Guardamos la ruta
for i in range(len(Hormigas)):
optimizer = AntColonyOptimizer(ants=Hormigas[i], evaporation_rate=.1, intensification=2, alpha=1, beta=1,
beta_evaporation_rate=0.5, choose_best=.1)
best = optimizer.fit(problem, 100) # Iniciamos el optimizador con 100 iteraciones
optimizer.plot()
MejorRuta.append(optimizer.best_path) #Agregamos la mejor ruta por cada numero de Hormigas
Beginning ACO Optimization with 100 iterations... Best score at iteration 0: 958526.71; overall: 958526.71 (0s) Best score at iteration 1: 833377.94; overall: 833377.94 (0s) Best score at iteration 2: 998090.98; overall: 833377.94 (0s) Best score at iteration 3: 1094889.29; overall: 833377.94 (0s) Best score at iteration 4: 1157027.73; overall: 833377.94 (0s) Best score at iteration 5: 911612.27; overall: 833377.94 (0s) Best score at iteration 6: 990091.12; overall: 833377.94 (0s) Best score at iteration 7: 1063159.76; overall: 833377.94 (0s) Best score at iteration 8: 878125.67; overall: 833377.94 (0s) Best score at iteration 9: 1105489.97; overall: 833377.94 (0s) Best score at iteration 10: 946840.45; overall: 833377.94 (0s) Best score at iteration 11: 850775.99; overall: 833377.94 (0s) Best score at iteration 12: 900991.92; overall: 833377.94 (0s) Best score at iteration 13: 1016693.63; overall: 833377.94 (0s) Best score at iteration 14: 824877.43; overall: 824877.43 (0s) Best score at iteration 15: 962244.56; overall: 824877.43 (0s) Best score at iteration 16: 865798.54; overall: 824877.43 (0s) Best score at iteration 17: 817484.4; overall: 817484.4 (0s) Best score at iteration 18: 755048.45; overall: 755048.45 (0s) Best score at iteration 19: 880523.13; overall: 755048.45 (0s) Best score at iteration 20: 867293.41; overall: 755048.45 (0s) Best score at iteration 21: 933027.84; overall: 755048.45 (0s) Best score at iteration 22: 874751.98; overall: 755048.45 (0s) Best score at iteration 23: 841491.38; overall: 755048.45 (0s) Best score at iteration 24: 783664.25; overall: 755048.45 (0s) Best score at iteration 25: 860620.6; overall: 755048.45 (0s) Best score at iteration 26: 736701.69; overall: 736701.69 (0s) Best score at iteration 27: 736701.69; overall: 736701.69 (0s) Best score at iteration 28: 782859.3; overall: 736701.69 (0s) Best score at iteration 29: 898906.44; overall: 736701.69 (0s) Best score at iteration 30: 778318.48; overall: 736701.69 (0s) Best score at iteration 31: 773751.46; overall: 736701.69 (0s) Best score at iteration 32: 747642.54; overall: 736701.69 (0s) Best score at iteration 33: 736701.69; overall: 736701.69 (0s) Best score at iteration 34: 736701.69; overall: 736701.69 (0s) Best score at iteration 35: 736701.69; overall: 736701.69 (0s) Best score at iteration 36: 736701.69; overall: 736701.69 (0s) Best score at iteration 37: 736701.69; overall: 736701.69 (0s) Best score at iteration 38: 736701.69; overall: 736701.69 (0s) Best score at iteration 39: 736701.69; overall: 736701.69 (0s) Best score at iteration 40: 736701.69; overall: 736701.69 (0s) Best score at iteration 41: 736701.69; overall: 736701.69 (0s) Best score at iteration 42: 736701.69; overall: 736701.69 (0s) Best score at iteration 43: 736701.69; overall: 736701.69 (0s) Best score at iteration 44: 736701.69; overall: 736701.69 (0s) Best score at iteration 45: 736701.69; overall: 736701.69 (0s) Best score at iteration 46: 736701.69; overall: 736701.69 (0s) Best score at iteration 47: 736701.69; overall: 736701.69 (0s) Best score at iteration 48: 736701.69; overall: 736701.69 (0s) Best score at iteration 49: 729205.86; overall: 729205.86 (0s) Best score at iteration 50: 729205.86; overall: 729205.86 (0s) Best score at iteration 51: 729205.86; overall: 729205.86 (0s) Best score at iteration 52: 729205.86; overall: 729205.86 (0s) Best score at iteration 53: 729205.86; overall: 729205.86 (0s) Best score at iteration 54: 729205.86; overall: 729205.86 (0s) Best score at iteration 55: 729205.86; overall: 729205.86 (0s) Best score at iteration 56: 729205.86; overall: 729205.86 (0s) Best score at iteration 57: 729205.86; overall: 729205.86 (0s) Best score at iteration 58: 729205.86; overall: 729205.86 (0s) Best score at iteration 59: 729205.86; overall: 729205.86 (0s) Best score at iteration 60: 729205.86; overall: 729205.86 (0s) Best score at iteration 61: 729205.86; overall: 729205.86 (0s) Best score at iteration 62: 729205.86; overall: 729205.86 (0s) Best score at iteration 63: 729205.86; overall: 729205.86 (0s) Best score at iteration 64: 729205.86; overall: 729205.86 (0s) Best score at iteration 65: 729205.86; overall: 729205.86 (0s) Best score at iteration 66: 729205.86; overall: 729205.86 (0s) Best score at iteration 67: 729205.86; overall: 729205.86 (0s) Best score at iteration 68: 729205.86; overall: 729205.86 (0s) Best score at iteration 69: 729205.86; overall: 729205.86 (0s) Best score at iteration 70: 729205.86; overall: 729205.86 (0s) Best score at iteration 71: 729205.86; overall: 729205.86 (0s) Best score at iteration 72: 729205.86; overall: 729205.86 (0s) Best score at iteration 73: 729205.86; overall: 729205.86 (0s) Best score at iteration 74: 729205.86; overall: 729205.86 (0s) Best score at iteration 75: 729205.86; overall: 729205.86 (0s) Best score at iteration 76: 729205.86; overall: 729205.86 (0s) Best score at iteration 77: 729205.86; overall: 729205.86 (0s) Best score at iteration 78: 729205.86; overall: 729205.86 (0s) Best score at iteration 79: 729205.86; overall: 729205.86 (0s) Best score at iteration 80: 729205.86; overall: 729205.86 (0s) Best score at iteration 81: 729205.86; overall: 729205.86 (0s) Best score at iteration 82: 729205.86; overall: 729205.86 (0s) Best score at iteration 83: 729205.86; overall: 729205.86 (0s) Best score at iteration 84: 729205.86; overall: 729205.86 (0s) Best score at iteration 85: 729205.86; overall: 729205.86 (0s) Best score at iteration 86: 729205.86; overall: 729205.86 (0s) Best score at iteration 87: 729205.86; overall: 729205.86 (0s) Best score at iteration 88: 729205.86; overall: 729205.86 (0s) Best score at iteration 89: 729205.86; overall: 729205.86 (0s) Best score at iteration 90: 729205.86; overall: 729205.86 (0s) Best score at iteration 91: 729205.86; overall: 729205.86 (0s) Best score at iteration 92: 729205.86; overall: 729205.86 (0s) Best score at iteration 93: 729205.86; overall: 729205.86 (0s) Best score at iteration 94: 729205.86; overall: 729205.86 (0s) Best score at iteration 95: 729205.86; overall: 729205.86 (0s) Best score at iteration 96: 729205.86; overall: 729205.86 (0s) Best score at iteration 97: 729205.86; overall: 729205.86 (0s) Best score at iteration 98: 729205.86; overall: 729205.86 (0s) Best score at iteration 99: 729205.86; overall: 729205.86 (0s) ACO fitted. Runtime: 0 minutes. Best score: 729205.8646025728
Beginning ACO Optimization with 100 iterations... Best score at iteration 0: 676258.8; overall: 676258.8 (0s) Best score at iteration 1: 1115739.0; overall: 676258.8 (0s) Best score at iteration 2: 848845.0; overall: 676258.8 (0s) Best score at iteration 3: 980641.06; overall: 676258.8 (0s) Best score at iteration 4: 929794.16; overall: 676258.8 (0s) Best score at iteration 5: 815519.3; overall: 676258.8 (0s) Best score at iteration 6: 937332.46; overall: 676258.8 (0s) Best score at iteration 7: 761930.4; overall: 676258.8 (0s) Best score at iteration 8: 833334.04; overall: 676258.8 (0s) Best score at iteration 9: 880173.73; overall: 676258.8 (0s) Best score at iteration 10: 737772.55; overall: 676258.8 (0s) Best score at iteration 11: 700355.48; overall: 676258.8 (0s) Best score at iteration 12: 806433.02; overall: 676258.8 (0s) Best score at iteration 13: 800018.72; overall: 676258.8 (0s) Best score at iteration 14: 853380.28; overall: 676258.8 (0s) Best score at iteration 15: 737772.55; overall: 676258.8 (0s) Best score at iteration 16: 767859.24; overall: 676258.8 (0s) Best score at iteration 17: 684849.0; overall: 676258.8 (0s) Best score at iteration 18: 720133.72; overall: 676258.8 (0s) Best score at iteration 19: 701458.38; overall: 676258.8 (0s) Best score at iteration 20: 684418.67; overall: 676258.8 (0s) Best score at iteration 21: 737772.55; overall: 676258.8 (0s) Best score at iteration 22: 709301.21; overall: 676258.8 (0s) Best score at iteration 23: 737772.55; overall: 676258.8 (0s) Best score at iteration 24: 735301.04; overall: 676258.8 (0s) Best score at iteration 25: 700355.48; overall: 676258.8 (0s) Best score at iteration 26: 685076.97; overall: 676258.8 (0s) Best score at iteration 27: 685076.97; overall: 676258.8 (0s) Best score at iteration 28: 708332.91; overall: 676258.8 (0s) Best score at iteration 29: 735883.45; overall: 676258.8 (0s) Best score at iteration 30: 684418.67; overall: 676258.8 (0s) Best score at iteration 31: 685076.97; overall: 676258.8 (0s) Best score at iteration 32: 676258.8; overall: 676258.8 (0s) Best score at iteration 33: 737772.55; overall: 676258.8 (0s) Best score at iteration 34: 684418.67; overall: 676258.8 (0s) Best score at iteration 35: 668281.37; overall: 668281.37 (0s) Best score at iteration 36: 668281.37; overall: 668281.37 (0s) Best score at iteration 37: 684418.67; overall: 668281.37 (0s) Best score at iteration 38: 684418.67; overall: 668281.37 (0s) Best score at iteration 39: 676258.8; overall: 668281.37 (0s) Best score at iteration 40: 668281.37; overall: 668281.37 (0s) Best score at iteration 41: 668281.37; overall: 668281.37 (0s) Best score at iteration 42: 668281.37; overall: 668281.37 (0s) Best score at iteration 43: 668281.37; overall: 668281.37 (0s) Best score at iteration 44: 668281.37; overall: 668281.37 (0s) Best score at iteration 45: 668281.37; overall: 668281.37 (0s) Best score at iteration 46: 668281.37; overall: 668281.37 (0s) Best score at iteration 47: 668281.37; overall: 668281.37 (0s) Best score at iteration 48: 668281.37; overall: 668281.37 (0s) Best score at iteration 49: 668281.37; overall: 668281.37 (0s) Best score at iteration 50: 668281.37; overall: 668281.37 (0s) Best score at iteration 51: 668281.37; overall: 668281.37 (0s) Best score at iteration 52: 668281.37; overall: 668281.37 (0s) Best score at iteration 53: 668281.37; overall: 668281.37 (0s) Best score at iteration 54: 668281.37; overall: 668281.37 (0s) Best score at iteration 55: 668281.37; overall: 668281.37 (0s) Best score at iteration 56: 668281.37; overall: 668281.37 (0s) Best score at iteration 57: 668281.37; overall: 668281.37 (0s) Best score at iteration 58: 668281.37; overall: 668281.37 (0s) Best score at iteration 59: 668281.37; overall: 668281.37 (0s) Best score at iteration 60: 668281.37; overall: 668281.37 (0s) Best score at iteration 61: 668281.37; overall: 668281.37 (0s) Best score at iteration 62: 668281.37; overall: 668281.37 (0s) Stopping early due to 20 iterations of the same score. ACO fitted. Runtime: 0 minutes. Best score: 668281.3692114921
Beginning ACO Optimization with 100 iterations... Best score at iteration 0: 725468.99; overall: 725468.99 (0s) Best score at iteration 1: 804441.6; overall: 725468.99 (0s) Best score at iteration 2: 824961.53; overall: 725468.99 (0s) Best score at iteration 3: 898382.85; overall: 725468.99 (0s) Best score at iteration 4: 897055.51; overall: 725468.99 (0s) Best score at iteration 5: 763540.39; overall: 725468.99 (0s) Best score at iteration 6: 754611.39; overall: 725468.99 (0s) Best score at iteration 7: 740694.62; overall: 725468.99 (0s) Best score at iteration 8: 739618.61; overall: 725468.99 (0s) Best score at iteration 9: 746949.15; overall: 725468.99 (0s) Best score at iteration 10: 740694.62; overall: 725468.99 (0s) Best score at iteration 11: 750253.08; overall: 725468.99 (0s) Best score at iteration 12: 727348.28; overall: 725468.99 (0s) Best score at iteration 13: 731547.84; overall: 725468.99 (0s) Best score at iteration 14: 727348.28; overall: 725468.99 (0s) Best score at iteration 15: 718515.62; overall: 718515.62 (0s) Best score at iteration 16: 727348.28; overall: 718515.62 (0s) Best score at iteration 17: 719899.96; overall: 718515.62 (0s) Best score at iteration 18: 684749.9; overall: 684749.9 (0s) Best score at iteration 19: 719671.99; overall: 684749.9 (0s) Best score at iteration 20: 684749.9; overall: 684749.9 (0s) Best score at iteration 21: 684749.9; overall: 684749.9 (0s) Best score at iteration 22: 684749.9; overall: 684749.9 (0s) Best score at iteration 23: 684749.9; overall: 684749.9 (0s) Best score at iteration 24: 668281.37; overall: 668281.37 (0s) Best score at iteration 25: 668281.37; overall: 668281.37 (0s) Best score at iteration 26: 668281.37; overall: 668281.37 (0s) Best score at iteration 27: 668281.37; overall: 668281.37 (0s) Best score at iteration 28: 668281.37; overall: 668281.37 (0s) Best score at iteration 29: 668281.37; overall: 668281.37 (0s) Best score at iteration 30: 668281.37; overall: 668281.37 (0s) Best score at iteration 31: 668281.37; overall: 668281.37 (0s) Best score at iteration 32: 668281.37; overall: 668281.37 (0s) Best score at iteration 33: 668281.37; overall: 668281.37 (0s) Best score at iteration 34: 668281.37; overall: 668281.37 (0s) Best score at iteration 35: 668281.37; overall: 668281.37 (0s) Best score at iteration 36: 668281.37; overall: 668281.37 (0s) Best score at iteration 37: 668281.37; overall: 668281.37 (0s) Best score at iteration 38: 668281.37; overall: 668281.37 (0s) Best score at iteration 39: 668281.37; overall: 668281.37 (0s) Best score at iteration 40: 668281.37; overall: 668281.37 (0s) Best score at iteration 41: 668281.37; overall: 668281.37 (0s) Best score at iteration 42: 668281.37; overall: 668281.37 (0s) Best score at iteration 43: 668281.37; overall: 668281.37 (0s) Best score at iteration 44: 668281.37; overall: 668281.37 (0s) Best score at iteration 45: 668281.37; overall: 668281.37 (0s) Best score at iteration 46: 668281.37; overall: 668281.37 (0s) Best score at iteration 47: 668281.37; overall: 668281.37 (0s) Best score at iteration 48: 668281.37; overall: 668281.37 (0s) Stopping early due to 20 iterations of the same score. ACO fitted. Runtime: 0 minutes. Best score: 668281.3692114921
Beginning ACO Optimization with 100 iterations... Best score at iteration 0: 775390.1; overall: 775390.1 (0s) Best score at iteration 1: 885962.49; overall: 775390.1 (0s) Best score at iteration 2: 726420.18; overall: 726420.18 (0s) Best score at iteration 3: 880384.95; overall: 726420.18 (0s) Best score at iteration 4: 799468.6; overall: 726420.18 (0s) Best score at iteration 5: 726420.18; overall: 726420.18 (0s) Best score at iteration 6: 766951.4; overall: 726420.18 (0s) Best score at iteration 7: 731997.72; overall: 726420.18 (0s) Best score at iteration 8: 726283.37; overall: 726283.37 (0s) Best score at iteration 9: 726420.18; overall: 726283.37 (0s) Best score at iteration 10: 724611.44; overall: 724611.44 (0s) Best score at iteration 11: 726420.18; overall: 724611.44 (0s) Best score at iteration 12: 726283.37; overall: 724611.44 (0s) Best score at iteration 13: 724611.44; overall: 724611.44 (0s) Best score at iteration 14: 715136.57; overall: 715136.57 (0s) Best score at iteration 15: 713699.13; overall: 713699.13 (0s) Best score at iteration 16: 715136.57; overall: 713699.13 (0s) Best score at iteration 17: 715136.57; overall: 713699.13 (0s) Best score at iteration 18: 715136.57; overall: 713699.13 (0s) Best score at iteration 19: 696726.59; overall: 696726.59 (0s) Best score at iteration 20: 715136.57; overall: 696726.59 (0s) Best score at iteration 21: 713587.74; overall: 696726.59 (0s) Best score at iteration 22: 696726.59; overall: 696726.59 (0s) Best score at iteration 23: 693705.55; overall: 693705.55 (0s) Best score at iteration 24: 675295.57; overall: 675295.57 (0s) Best score at iteration 25: 666529.07; overall: 666529.07 (0s) Best score at iteration 26: 684939.04; overall: 666529.07 (0s) Best score at iteration 27: 666529.07; overall: 666529.07 (0s) Best score at iteration 28: 666529.07; overall: 666529.07 (0s) Best score at iteration 29: 666529.07; overall: 666529.07 (0s) Best score at iteration 30: 666529.07; overall: 666529.07 (0s) Best score at iteration 31: 666529.07; overall: 666529.07 (0s) Best score at iteration 32: 666529.07; overall: 666529.07 (0s) Best score at iteration 33: 666392.26; overall: 666392.26 (0s) Best score at iteration 34: 666392.26; overall: 666392.26 (0s) Best score at iteration 35: 666529.07; overall: 666392.26 (0s) Best score at iteration 36: 666392.26; overall: 666392.26 (0s) Best score at iteration 37: 666392.26; overall: 666392.26 (0s) Best score at iteration 38: 666392.26; overall: 666392.26 (0s) Best score at iteration 39: 666392.26; overall: 666392.26 (0s) Best score at iteration 40: 666392.26; overall: 666392.26 (0s) Best score at iteration 41: 666392.26; overall: 666392.26 (0s) Best score at iteration 42: 666392.26; overall: 666392.26 (0s) Best score at iteration 43: 666392.26; overall: 666392.26 (0s) Best score at iteration 44: 666392.26; overall: 666392.26 (0s) Best score at iteration 45: 666392.26; overall: 666392.26 (0s) Best score at iteration 46: 666392.26; overall: 666392.26 (0s) Best score at iteration 47: 666392.26; overall: 666392.26 (0s) Best score at iteration 48: 666392.26; overall: 666392.26 (0s) Best score at iteration 49: 666392.26; overall: 666392.26 (0s) Best score at iteration 50: 666392.26; overall: 666392.26 (0s) Best score at iteration 51: 666392.26; overall: 666392.26 (0s) Best score at iteration 52: 666392.26; overall: 666392.26 (0s) Best score at iteration 53: 666392.26; overall: 666392.26 (0s) Best score at iteration 54: 666392.26; overall: 666392.26 (0s) Best score at iteration 55: 666392.26; overall: 666392.26 (0s) Stopping early due to 20 iterations of the same score. ACO fitted. Runtime: 0 minutes. Best score: 666392.259629674
Beginning ACO Optimization with 100 iterations... Best score at iteration 0: 857230.94; overall: 857230.94 (0s) Best score at iteration 1: 865477.55; overall: 857230.94 (0s) Best score at iteration 2: 888071.88; overall: 857230.94 (0s) Best score at iteration 3: 916103.29; overall: 857230.94 (0s) Best score at iteration 4: 825349.44; overall: 825349.44 (0s) Best score at iteration 5: 915775.09; overall: 825349.44 (0s) Best score at iteration 6: 840176.01; overall: 825349.44 (0s) Best score at iteration 7: 833182.93; overall: 825349.44 (0s) Best score at iteration 8: 687450.69; overall: 687450.69 (0s) Best score at iteration 9: 765471.01; overall: 687450.69 (0s) Best score at iteration 10: 695595.78; overall: 687450.69 (0s) Best score at iteration 11: 709547.32; overall: 687450.69 (0s) Best score at iteration 12: 756948.6; overall: 687450.69 (0s) Best score at iteration 13: 695342.3; overall: 687450.69 (0s) Best score at iteration 14: 693615.5; overall: 687450.69 (0s) Best score at iteration 15: 730712.44; overall: 687450.69 (0s) Best score at iteration 16: 703433.44; overall: 687450.69 (0s) Best score at iteration 17: 693478.7; overall: 687450.69 (0s) Best score at iteration 18: 693615.5; overall: 687450.69 (0s) Best score at iteration 19: 688284.18; overall: 687450.69 (0s) Best score at iteration 20: 693615.5; overall: 687450.69 (0s) Best score at iteration 21: 679161.08; overall: 679161.08 (0s) Best score at iteration 22: 680124.51; overall: 679161.08 (0s) Best score at iteration 23: 688804.55; overall: 679161.08 (0s) Best score at iteration 24: 679161.08; overall: 679161.08 (0s) Best score at iteration 25: 670394.58; overall: 670394.58 (0s) Best score at iteration 26: 670394.58; overall: 670394.58 (0s) Best score at iteration 27: 670257.77; overall: 670257.77 (0s) Best score at iteration 28: 670257.77; overall: 670257.77 (0s) Best score at iteration 29: 670257.77; overall: 670257.77 (0s) Best score at iteration 30: 670257.77; overall: 670257.77 (0s) Best score at iteration 31: 670257.77; overall: 670257.77 (0s) Best score at iteration 32: 670257.77; overall: 670257.77 (0s) Best score at iteration 33: 670257.77; overall: 670257.77 (0s) Best score at iteration 34: 670257.77; overall: 670257.77 (0s) Best score at iteration 35: 670257.77; overall: 670257.77 (0s) Best score at iteration 36: 670257.77; overall: 670257.77 (0s) Best score at iteration 37: 670257.77; overall: 670257.77 (0s) Best score at iteration 38: 670257.77; overall: 670257.77 (0s) Best score at iteration 39: 670257.77; overall: 670257.77 (0s) Best score at iteration 40: 670257.77; overall: 670257.77 (0s) Best score at iteration 41: 670257.77; overall: 670257.77 (0s) Best score at iteration 42: 670257.77; overall: 670257.77 (0s) Best score at iteration 43: 670257.77; overall: 670257.77 (0s) Best score at iteration 44: 670257.77; overall: 670257.77 (0s) Best score at iteration 45: 670257.77; overall: 670257.77 (0s) Best score at iteration 46: 670257.77; overall: 670257.77 (0s) Best score at iteration 47: 670257.77; overall: 670257.77 (0s) Best score at iteration 48: 670257.77; overall: 670257.77 (0s) Best score at iteration 49: 670257.77; overall: 670257.77 (0s) Best score at iteration 50: 670257.77; overall: 670257.77 (0s) Best score at iteration 51: 670257.77; overall: 670257.77 (0s) Stopping early due to 20 iterations of the same score. ACO fitted. Runtime: 0 minutes. Best score: 670257.7702044595
Calculamos el costo total del recorrido
def costo_total(ruta, costos):
costo = 0
for i in range(len(ruta) - 1):
costo += MCTC[ruta[i]][ruta[i+1]]
return costo
for i in range(len(MejorRuta)):
print("El costo total de la ruta es: $", round(costo_total(MejorRuta[i], MCTC), 2), " COP para un total de ", Hormigas[i], "Hormigas")
El costo total de la ruta es: $ 729205.86 COP para un total de 10 Hormigas El costo total de la ruta es: $ 668281.37 COP para un total de 20 Hormigas El costo total de la ruta es: $ 668281.37 COP para un total de 80 Hormigas El costo total de la ruta es: $ 666392.26 COP para un total de 100 Hormigas El costo total de la ruta es: $ 670257.77 COP para un total de 120 Hormigas
#Cada vez que ejecutemos esta celda y el script en general borramos los archivos para que no se acumulen plot de otras seseiones
folder_path = "GifVCH"
for file_name in os.listdir(folder_path):
file_path = os.path.join(folder_path, file_name)
if os.path.isfile(file_path):
os.remove(file_path)
# Lista para almacenar las rutas de las rutas de las imágenes PNG
ruta_imagenes = []
for i in range(len(MejorRuta)):
texto = "La ruta es: "+"\n" # texto inical para añadir las rutas al plot
texto2 = "El costo total de la ruta es: $" + str(round(costo_total(MejorRuta[i], MCTC), 2))+"COP" # Texto para añadir el valor de la ruta
current_sol = MejorRuta[i] # Obtenbemos la soulcion por iteracion
for j in range(len(current_sol)):
texto += str(ciudades[current_sol[j]])+"\n"
# Aqui obtenemos la ruta segun sea la solucion
rutacol = [df.iloc[i]['geometry']for i in MejorRuta[i]]
#creamos el plot de la ruta
fig, ax = plt.subplots(figsize=(14, 14))
mapa.plot(ax=ax, color='lime', edgecolor='black')
df.plot(ax=ax, markersize=50, color='red')
for i, txt in enumerate(df['Ciudad']):
ax.annotate(txt, (df['geometry'][i].x,
df['geometry'][i].y), fontsize=12)
df.plot(ax=plt.gca(), color='red', markersize=50, zorder=3)
plt.plot([p.x for p in rutacol], [p.y for p in rutacol],
color='blue', linestyle='--', linewidth=2, dashes=[2, 2], zorder=2)
plt.text(250000, 1000000, texto, ha='center', va='center', size=14,
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.9))
plt.text(550000, 2000000, texto2, ha='center', va='center', size=14,
bbox=dict(boxstyle='round', facecolor='cyan', alpha=0.9))
plt.title("Problema del viajero: Colonias de Hormigas")
# Generar un nombre de archivo único con un número aleatorio
nombre_archivo = f"GifVCH/imagen_{i}_{random.randint(0, 100000)}.png"
texto = ""
# Agregar la ruta de la imagen PNG a la lista
ruta_imagenes.append(nombre_archivo)
# Guardar la figura como imagen PNG
plt.savefig(nombre_archivo)
# Crear archivo GIF a partir de las imágenes PNG
with imageio.get_writer('GifVCH/animacionCH.gif', mode='I', duration=0.8) as writer:
for ruta in ruta_imagenes:
image = imageio.imread(ruta)
writer.append_data(image)
#Variables necesarias para la creacion de gif
num_imagenes = 10
imagenes = []
def tsp_fitness(solution, s_idx):
#Realizamos unos ajustes a soltion puesto que no se encuentran en formato que nos deje iterar
solution = solution.tolist()
solution = [int(x) for x in solution]
# Verificar que no haya ciudades repetidas en la solución
if len(set(solution)) != len(solution):
return 0 # retorna aptitud cero si hay ciudades repetidas en la solución
# Verificar que la solución contenga todas las ciudades
if set(solution) != set(range(len(ciudades))):
return 0 # retorna aptitud cero si falta alguna ciudad en la solución
# Calcular el costo total del recorrido
distance = 0
for i in range(len(solution)-1):
distance += MCTC[solution[i]][solution[i+1]]
distance += MCTC[solution[-1]][solution[0]] # vuevle al punto de inicio
# Retornar la inversa de la distancia como aptitud (para maximizar)
return 1/distance
Generaciones = [100, 150, 270, 500, 1120, 1500]
mejor_soluciones = []
sol_per_pop = 15
soluciones_iniciales = [list(set(np.random.permutation(idx_ciudades)))
for _ in range(sol_per_pop)]
for i in range(len(Generaciones)):
num_generations = Generaciones[i]
# Definir el modelo de optimización genética con PyGAD
ViajeroGA = pygad.GA(num_generations=num_generations,
num_parents_mating=4,
sol_per_pop=sol_per_pop,
num_genes=len(ciudades),
initial_population=soluciones_iniciales,
mutation_percent_genes=30,
fitness_func=tsp_fitness,
mutation_type="scramble",
init_range_low=1,
init_range_high=15,
parent_selection_type="rank")
# Ejecutar la optimización genética
ViajeroGA.run()
# Obtener la mejor solución encontrada
solution, solution_fitness, solution_idx = ViajeroGA.best_solution()
solution_indices = [int(x) for x in solution.tolist()]
mejor_soluciones.append(solution_indices)
MejorCiudades = []
for i in solution_indices:
MejorCiudades.append(ciudades[i])
# Imprimir la mejor solución encontrada
print("Mejor solución encontrada en:",
num_generations, "iteraciones es:", MejorCiudades)
Mejor solución encontrada en: 100 iteraciones es: ['Cúcuta', 'Bucaramanga', 'Bogota', 'Armenia', 'Tuluá', 'Pasto', 'Palmira', 'Pereira', 'Manizales', 'Medellín', 'Montería', 'Cartagena', 'Soledad', 'Barranquilla', 'Valledupar'] Mejor solución encontrada en: 150 iteraciones es: ['Montería', 'Bogota', 'Manizales', 'Medellín', 'Tuluá', 'Palmira', 'Pasto', 'Pereira', 'Armenia', 'Valledupar', 'Barranquilla', 'Soledad', 'Cartagena', 'Bucaramanga', 'Cúcuta'] Mejor solución encontrada en: 270 iteraciones es: ['Manizales', 'Armenia', 'Pereira', 'Tuluá', 'Palmira', 'Pasto', 'Bogota', 'Bucaramanga', 'Cúcuta', 'Valledupar', 'Cartagena', 'Barranquilla', 'Soledad', 'Montería', 'Medellín'] Mejor solución encontrada en: 500 iteraciones es: ['Bogota', 'Armenia', 'Tuluá', 'Palmira', 'Pasto', 'Pereira', 'Manizales', 'Medellín', 'Montería', 'Cúcuta', 'Cartagena', 'Soledad', 'Barranquilla', 'Valledupar', 'Bucaramanga'] Mejor solución encontrada en: 1120 iteraciones es: ['Bogota', 'Armenia', 'Palmira', 'Pasto', 'Tuluá', 'Pereira', 'Medellín', 'Montería', 'Soledad', 'Barranquilla', 'Cartagena', 'Valledupar', 'Cúcuta', 'Bucaramanga', 'Manizales'] Mejor solución encontrada en: 1500 iteraciones es: ['Valledupar', 'Bogota', 'Armenia', 'Pasto', 'Palmira', 'Tuluá', 'Pereira', 'Manizales', 'Medellín', 'Montería', 'Barranquilla', 'Soledad', 'Cartagena', 'Cúcuta', 'Bucaramanga']
Calculamos el costo total del recorrido
def costo_total(ruta, costos):
costo = 0
for i in range(len(ruta) - 1):
costo += MCTC[ruta[i]][ruta[i+1]]
return costo
for i in range(len(mejor_soluciones)):
print("El costo total de la ruta es: $",round(costo_total(mejor_soluciones[i],MCTC),2)," COP para un total de ", Generaciones[i], "Iteraciones")
El costo total de la ruta es: $ 648272.0 COP para un total de 100 Iteraciones El costo total de la ruta es: $ 749751.41 COP para un total de 150 Iteraciones El costo total de la ruta es: $ 653715.53 COP para un total de 270 Iteraciones El costo total de la ruta es: $ 655548.19 COP para un total de 500 Iteraciones El costo total de la ruta es: $ 641403.75 COP para un total de 1120 Iteraciones El costo total de la ruta es: $ 662283.01 COP para un total de 1500 Iteraciones
#Cada vez que ejecutemos esta celda y el script en general borramos los archivos para que no se acumulen plot de otras seseiones
folder_path = "GifVGA"
for file_name in os.listdir(folder_path):
file_path = os.path.join(folder_path, file_name)
if os.path.isfile(file_path):
os.remove(file_path)
# Lista para almacenar las rutas de las rutas de las imágenes PNG
ruta_imagenes = []
for i in range(len(mejor_soluciones)):
texto = "La ruta es: "+"\n" # texto inical para añadir las rutas al plot
texto2 = "El costo total de la ruta es: $"+str(round(costo_total(mejor_soluciones[i], MCTC), 2))+"COP" #Texto para añadir el valor de la ruta
current_sol = mejor_soluciones[i] #Obtenbemos la soulcion por iteracion
for j in range(len(current_sol)):
texto += str(ciudades[current_sol[j]])+"\n"
rutacol = [df.iloc[i]['geometry']for i in mejor_soluciones[i]] #Aqui obtenemos la ruta segun sea la solucion
#creamos el plot de la ruta
fig, ax = plt.subplots(figsize=(14, 14))
mapa.plot(ax=ax, color='lime', edgecolor='black')
df.plot(ax=ax, markersize=50, color='red')
for i, txt in enumerate(df['Ciudad']):
ax.annotate(txt, (df['geometry'][i].x,
df['geometry'][i].y), fontsize=12)
df.plot(ax=plt.gca(), color='red', markersize=50, zorder=3)
plt.plot([p.x for p in rutacol], [p.y for p in rutacol],
color='blue', linestyle='--', linewidth=2, dashes=[2, 2], zorder=2)
plt.text(250000, 1000000, texto, ha='center', va='center', size=14,
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.9))
plt.text(550000, 2000000, texto2, ha='center', va='center', size=14,
bbox=dict(boxstyle='round', facecolor='cyan', alpha=0.9))
plt.title("Problema del viajero: Colonias de Hormigas")
# Generar un nombre de archivo único con un número aleatorio
nombre_archivo = f"GifVGA/imagen_{i}_{random.randint(0, 100000)}.png"
texto=""
# Agregar la ruta de la imagen PNG a la lista
ruta_imagenes.append(nombre_archivo)
# Guardar la figura como imagen PNG
plt.savefig(nombre_archivo)
# Crear archivo GIF a partir de las imágenes PNG
with imageio.get_writer('GifVGA/animacionGA.gif', mode='I', duration=0.8) as writer:
for ruta in ruta_imagenes:
image = imageio.imread(ruta)
writer.append_data(image)
from IPython.display import Image, display
gif1 = Image(filename='GifVCH/animacionCH.gif')
gif2 = Image(filename='GifVGA/animacionGA.gif')
display(gif1, gif2)
<IPython.core.display.Image object>
<IPython.core.display.Image object>
La optimizacion del problema por Colonias de hormigas nos da las rutas mas baratas, tambien se nota que por este metodo convergen mas rapido las soluciones a rutas Optimas